Matrix multiplication to get covariance matrix

10 visualizzazioni (ultimi 30 giorni)
Tintin Milou
Tintin Milou il 12 Mag 2021
Commentato: James Tursa il 18 Mag 2021
Hi,
Is there any chance that the following code could be optimized? The matrices are not sparse. I wonder whether the fact that the matrix A doesn't change might make it possible to pre-calculate an expression before the loop to avoid the double multiplication in this A*V*A' term.
Notes :
Since B is symmetric, V will be symmetric, too.
The matrix A can be rank-deficient and therefore V will also be rank-deficient.
Thank you!
T = 500;
N = 1200;
Vt = zeros(N,N,T);
B = randn(N,N);
B = (B+B')/2;
A = randn(N,N);
for t = 1:T
Vt(:,:,t+1) = A*Vt(:,:,t)*A'+B;
end

Risposte (2)

Jan
Jan il 12 Mag 2021
No, there is no chance for a significant speedup or vectorization, because the result of the former iteration is the input of the current one. One tiny improvement is to pre-allocate Vt to the exact matching size:
Vt = zeros(N, N, T + 1);
% ^^^
  4 Commenti
Walter Roberson
Walter Roberson il 14 Mag 2021
format long g
A = randn(5)
A = 5×5
1.63008317959338 -0.917913331834656 -0.240134840535616 -1.41517193395925 -0.55948141889238 -1.11684857867162 -0.316475244146763 -0.202510039340634 1.00697516701369 0.27010155043662 0.724379290944823 -0.253694739742835 1.73238399447882 -0.671399324869163 -1.47295169851759 -2.59262795755675 0.150174708929094 1.18753962293291 -2.46296756702745 0.447970415690913 1.61893080922682 -0.32335613338683 -1.53128636093959 1.16901375411328 0.0685926510984651
B = randn(5); B = (B+B')/2 %symmetric
B = 5×5
0.256366520738895 -0.264727426083824 -1.77963279150492 0.460091526428723 1.00166819378983 -0.264727426083824 -0.630122678243637 0.296024585465261 0.0920082912890164 0.185525817419845 -1.77963279150492 0.296024585465261 -0.183550726829703 -0.1631387939754 0.137857032974269 0.460091526428723 0.0920082912890164 -0.1631387939754 0.979340158549388 -0.602329849032058 1.00166819378983 0.185525817419845 0.137857032974269 -0.602329849032058 -1.39020920985583
Vt = B;
Vt = A*Vt*A' + B
Vt = 5×5
-0.309619613568741 1.0052304350345 -11.2579060646595 -0.558371010417369 5.21606220918395 1.0052304350345 -2.44932037312056 6.93199118768203 -1.74284563418364 -1.10370885072686 -11.2579060646595 6.93199118768203 -11.7750226767206 10.5900263787199 -4.12908870717582 -0.55837101041737 -1.74284563418364 10.5900263787199 25.2977711902404 -18.8861734546075 5.21606220918395 -1.10370885072686 -4.12908870717582 -18.8861734546075 11.8558329547431
Vt - Vt'
ans = 5×5
0 0 0 1.11022302462516e-15 -8.88178419700125e-16 0 0 2.66453525910038e-15 8.88178419700125e-16 -6.66133814775094e-16 0 -2.66453525910038e-15 0 -3.5527136788005e-15 0 -1.11022302462516e-15 -8.88178419700125e-16 3.5527136788005e-15 0 -7.105427357601e-15 8.88178419700125e-16 6.66133814775094e-16 0 7.105427357601e-15 0
Not symmetric.
Tintin Milou
Tintin Milou il 14 Mag 2021
Ok, Walter. I think theoretically, Vt is symmetric. It's just that Matlab's calculations are not 100% precise. For my purpose, I do not care about any differences from zero of that magnitude. I noticed that transforming the matrices to 'single' helps quite a bit in terms of speed and it is precise enough for my purpose. But that's a different issue.

Accedi per commentare.


James Tursa
James Tursa il 17 Mag 2021
Modificato: James Tursa il 17 Mag 2021
@Tintin Milou MATLAB variables do not have a "symmetric" attribute flag that I know of. Just because you create a variable that is symmetic with an expression such as (A+A') doesn't mean that MATLAB will necessarily know downstream that the matrix is symmetric. To MATLAB this result is simply just another matrix, and any subsequent operation downstream that depended on symmetry would have to check all the elements. Some functions will do this check but others will not.
There are certain matrix multiplication expressions that MATLAB can recognize where the result will necessarily be symmetric, and in those cases MATLAB can call symmetric BLAS routines in the background to speed up calculations and obtain an exact symmetric result. E.g.,
A = an arbitrary real matrix
C = A' * A is recognized by MATLAB as being symmetric and it will call a symmetric BLAS routine in the background.
The C result will take less time and the result is guaranteed to be exactly symmetric.
B = A'
D = B * A is not recognized by MATLAB as being symmetric, so a generic BLAS routine will be used.
The D result will take more time and the result is not guaranteed to be exactly symmetric. MATLAB doesn't keep track of where B came from and doesn't know it is exactly the transpose of A, so it doesn't do any symmetric optimizations for the operation.
And for the different expression A * B * A' where B is symmetric, MATLAB does not check that B is symmetric so it doesn't know that the result will be symmetric. There are no symmetric optimizations taken for this case, and generic matrix multiplies are done instead with the result not guaranteed to be exactly symmetric because of floating point arithmetic effects.
  2 Commenti
Jan
Jan il 18 Mag 2021
Modificato: Jan il 18 Mag 2021
It would be nice to have a set of flags to mark a matrix as symmetric or transposed. Because MATLAB is the "Matrix Laboratory" and the underlying BLAS and LAPACK functions support at least the transposed property and the implementation is cheap, I'd actually expected such flags. It looks like the JIT can omit an explicit transposition on the fly also.
James Tursa
James Tursa il 18 Mag 2021
There does appear to be several unused bits in the mxArray flags int field that could be used for this. E.g., an isSymmetric and/or isHermitian bit. Functions could examine this and force exact symmetric results when appropriate. The main complication might be the parser recognizing expressions that result in symmetry and appropriately setting the bit(s).

Accedi per commentare.

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by