Multiplication of submatrices of a matrix

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!

[Paul]

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.

Answer 1

[David Goodmanson]

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.

[Laura Hurt]

Hi @David Goodmanson. Thank you so much for answering my question. You exactly got the idea. Anewproduct is indeed the outcome I am looking for. However, l was hoping to get rid of a loop as in my case n is around 170 000. If I can obtain Anewblock with matrix operations I think it'd be great.

Below I show how I am tackling this. In particular, the line to create A_new uses repmat(c,n,1), which it is taking a long time to compute. I would like to take advantage of the upper triangular matrix feature of c to make the script much faster.

%This is what I am currently doing (borrowing your code):

j = 3;   % for example
n = 5;
k = 7;
c = triu(2*rand(j,j)-1);
b = 2*rand(j,k-j,n)-1;     % stack the b_i in the third dimension

%computationally expensive target:
for m=1:n
    A_new_target(j*(m - 1) + 1:j*m,:) = [c c*b(:,:,m) ];
end
A_new_product_target = A_new_target'*A_new_target;


%Less computationally intensive alternative
b_horizontal = reshape(permute(b, [1, 3, 2]),j ,[] );
cb_horizontal = c*b_horizontal;
cb_horizontal_3d = reshape(cb_horizontal, j, k-j, n);
A_new_temp = reshape(permute(cb_horizontal, [1,3,2]), n*j, k-j)
A_new = [repmat(c,n,1) A_new_temp];
Anewproduct = A_new'*A_new


%comparison
max(max(A_new_product_target - Anewproduct))

[David Goodmanson]

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.

Answer 2

[Sudarsanan A K]

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:

1. Upper Left Block: This is the sum of with itself N times, which is .
2. Upper Right Block: This involves the sum of products of with each , which simplifies to .
3. Lower Left Block: This is the transpose of the upper right block, so it's .
4. 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 = ');

A_new' * A_new = 

disp(A_newT_A_new);

    1.4028    1.6114    0.6523    2.1441    2.2367
    1.6114    1.8535    0.8302    2.5049    2.6392
    0.6523    0.8302    3.0882    2.4391    3.4435
    2.1441    2.5049    2.4391    4.5182    4.4340
    2.2367    2.6392    3.4435    4.4340    5.8354

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:');

Direct computation of A_new' * A_new:

disp(A_newT_A_new_direct);

    1.4028    1.6114    0.6523    2.1441    2.2367
    1.6114    1.8535    0.8302    2.5049    2.6392
    0.6523    0.8302    3.0882    2.4391    3.4435
    2.1441    2.5049    2.4391    4.5182    4.4340
    2.2367    2.6392    3.4435    4.4340    5.8354

% 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

The results are the same within the specified tolerance.

I hope this helps!

[Laura Hurt]

Thank you @Sudarsanan A K for such detailed answered. As I mentioned in my previous comment to @David Goodmanson's answer, I am struggling with the loops involved in the operations as the matrices dimensions that i am working with are large. For instance, N is around 170k, J is 146 and K = 176. For this reason the loop to compute the lower right block of in particular is too costly for me.