Multiplication of submatrices of a matrix
3 views (last 30 days)
Show older comments
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 Comment
Paul
on 13 Jun 2024
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.
Accepted Answer
David Goodmanson
on 14 Jun 2024
Edited: David Goodmanson
on 14 Jun 2024
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 Comments
David Goodmanson
on 14 Jun 2024
Edited: David Goodmanson
on 14 Jun 2024
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.
More Answers (1)
Sudarsanan A K
on 14 Jun 2024
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.
Given:
Task: Express in a simplified form.
First, note that When we compute , we are essentially doing a block matrix multiplication. Let us break it down:
- 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:
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.
% 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:
% 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!
See Also
Categories
Find more on Matrices and Arrays in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!