Multiplication of submatrices of a matrix
4 次查看(过去 30 天)
显示 更早的评论
I have a J X J matrix $C$ that is upper triangular. Also, C'C is positive definite. I also have a matrix $A$ formed by submatrices of size J X K as follows, A = [A_1 ; A_2 ; A_3 ; ... A_N ]
where A_i = [I B_i] for each $i$. I is the identity matrix of size J and B_i is a matrix of size J X (K-J). I create a new matrix A_new such that
A_new = [C*A_1; C*A_2,...;C*A_N]
Notice that A_new can be written as A_new = [C CB_1 ; C CB_2 ; ... C CB_N ].
Is there a way for me to write A_new' * A_new in terms of matrix operations or concatenations?
I have a loop in Matlab where I need to obtain a matrix A_new for different values of C. This is extremely costly computationally as the matrices B_i and C are huge. Any help is appreciated!
Thank you!
1 个评论
Paul
2024-6-13
Hi Laura,
I suggest uploading a small example of C and A in a .mat file (use the paperclip icon on the insert menu) and edit the question to include the code to construct A_new.
采纳的回答
David Goodmanson
2024-6-14
编辑:David Goodmanson
2024-6-14
Hi Laura, if I interpret this correctly, you have n rows (I'll use small letters mostly) of the form
Anew =
[c c*b1
c c*b2
.
.
c c*bn]
so that A has dimensions j*n x k, and the product Anew' * Anew will be k x k. One way to reduce the number of operations is
j = 3; % for example
n = 5;
k = 7;
c = triu(2*rand(j,j)-1);
b = 2*rand(j,k-j,n)-1; % the b_i are stacked in the third dimension
% reduced calculations
cc = c'*c;
bsum = sum(b,3);
Anewblock = b(:,:,1)'*cc*b(:,:,1);
for m = 2:n
Anewblock = Anewblock + b(:,:,m)'*cc*b(:,:,m);
end
Anewproduct = [n*cc cc*bsum % 2x2 block of matrices
bsum'*cc Anewblock]
% check using Anew
Anew = zeros(j*n,k);
for m = 1:n
Anew((1:j)+(m-1)*j,:) = [c c*b(:,:,m)];
end
Anewproductcheck = Anew'*Anew
max(max(Anewproduct -Anewproductcheck )) % should be small
For convenience here the b_i matrices are saved in the third dimension. You may or may not have enough memory to do that, but one way or another you have to sum the b_i matrices.
The lower right block of Anewproduct is not the most convenient thing, a sum of b_i * cc * b_i.
The upper right and lower left blocks are transposes of each other with ' and you may want to use that fact.
2 个评论
David Goodmanson
2024-6-14
编辑:David Goodmanson
2024-6-14
Hi Laura,
170k is a lot, but it isn't always the case that matrix operations are way faster. Some matrix operations are slow, such as transpose I believe, repmat it looks like from your experience, and others. At least one is fast (but I can't remember which one!). Anyway, it could happen that going to matrix operations does not lead to a large speed increase.
更多回答(1 个)
Sudarsanan A K
2024-6-14
Hi Laura,
Given the structure of your matrices, we can indeed express
in terms of matrix operations that might help you avoid direct computation through loops, especially for large matrices. Let us break down the components based on the information you have provided.
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1715476/image.png)
Given:
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1715481/image.png)
Task: Express
in a simplified form.
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1715486/image.png)
First, note that
When we compute
, we are essentially doing a block matrix multiplication. Let us break it down:
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1715491/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1715496/image.png)
- Upper Left Block: This is the sum of
with itself N times, which is
.
- Upper Right Block: This involves the sum of products of
with each
, which simplifies to
.
- Lower Left Block: This is the transpose of the upper right block, so it's
.
- Lower Right Block: This is a bit more complex, involving sums of products of
with each
. This results in a sum of matrices of the form
for all
. However, for simplification, we can consider it as the sum of each
for all i, assuming the cross terms are not directly simplified without additional properties or context.
Given these components,
can be expressed as:
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1715556/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1715561/image.png)
Here is a simple MATLAB demonstration:
We will assume you have the matrices
stored in a three-dimensional array (where each "slice" of the array along the third dimension is one of the
matrices), and you are computing
for a specific C matrix.
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1715566/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1715571/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1715576/image.png)
% Example dimensions
J = 3; % Dimension of C and I
K = 5; % Dimension of A_i
N = 4; % Number of A_i matrices
% Example C matrix (upper triangular and ensuring C'C is positive definite)
C = triu(rand(J, J));
while det(C'*C) <= 0
C = triu(rand(J, J)); % Ensure C'C is positive definite
end
% Generating example B_i matrices
B = rand(J, K-J, N); % 3D array to hold B_i matrices
% Precompute sums involving B_i for efficiency
sumBi = zeros(J, K-J);
for i = 1:N
sumBi = sumBi + B(:, :, i);
end
% Compute components for A_new' * A_new
CC = C' * C;
sumBiPrimeCC = sumBi' * CC;
CCsumBi = CC * sumBi;
% Initialize the lower right block
lowerRightBlock = zeros(K-J, K-J);
for i = 1:N
lowerRightBlock = lowerRightBlock + (C * B(:, :, i))' * (C * B(:, :, i));
end
% Assemble A_new' * A_new
A_newT_A_new = [N * CC, CCsumBi; sumBiPrimeCC, lowerRightBlock];
% Display the result
disp('A_new'' * A_new = ');
disp(A_newT_A_new);
To validate the computation of
using the proposed method, you can compare it against a direct computation of
where
Here is how you can do it in MATLAB:
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1715581/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1715586/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1715591/image.png)
% Direct computation of A_new
A_new_direct = [];
for i = 1:N
A_i = [eye(J), B(:, :, i)];
A_new_direct = [A_new_direct; C * A_i];
end
A_newT_A_new_direct = A_new_direct' * A_new_direct;
% Display the direct computation result
disp('Direct computation of A_new'' * A_new:');
disp(A_newT_A_new_direct);
% Assuming A_newT_A_new was computed using the proposed method (previous code)
% Comparison
difference = A_newT_A_new - A_newT_A_new_direct;
tolerance = 1e-10; % Tolerance for numerical comparison
if max(abs(difference(:))) < tolerance
disp('The results are the same within the specified tolerance.');
else
disp('The results differ.');
end
I hope this helps!
另请参阅
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!