Symmetric Matrix-Vector Multiplication with only lower triangular stored
32 次查看(过去 30 天)
显示 更早的评论
I have a very big (n>1000000) sparse symmetric positive definite matrix A and I want to multiply it efficiently by a vector x. Only the lower triangular of matrix A is stored with function "sparse". Is there any function inbuilt or external that anyone knows of that takes as input only the lower triangular and performs the complete operation Ax without having to recover whole A (adding the strict upper triangular again)?
Thank you very much,
Jan
0 个评论
采纳的回答
Jan
2021-4-24
编辑:Jan
2021-4-24
Plese check if this is efficient in your case:
A = rand(5, 5);
A = A + A'; % Symmetric example matrix
B = tril(A); % Left triangular part
x = rand(5, 1);
y1 = A * x;
y2 = B * x + B.' * x - diag(B) .* x;
y1 - y2
I assume that Matlab's JIT can handle B.' * x without computing B.' explicitly. Alternative:
y2 = B * x + B.' * x - diag(diag(B)) * x;
6 个评论
Clayton Gotberg
2021-4-25
You may try taking advantage of the reshape command as @Bruno Luong suggested below, to see if it's any faster than the transpose for row and column vectors.
When you say that the difference between A*x and B*x+B'x-diag(B).*x gets larger, do you mean that the time cost is larger or that the numbers actually change? Mathematically, that shouldn't be happening, but computers+numbers == funny things.
If the time cost just becomes too high, I can take another pass at it with an eye to exploiting the symmetry of the problem. My gut reaction was this handy identity, but there have to be a few dozen more ways to skin this cat.
更多回答(1 个)
Clayton Gotberg
2021-4-24
If
is a symmetric matrix and
is the lower triangular part of the matrix and
is the upper triangular part of the matrix:
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/595600/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/595605/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/595610/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/595615/image.png)
where the diagonal function only finds the diagonal elements of
. This is because of a few relations:
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/595620/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/595625/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/595630/image.png)
To save time and space on MATLAB (because the upper triangular matrix will take up much more space), take advantage of the relations:
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/595635/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/595640/image.png)
To get:
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/595645/image.png)
Now, the MATLAB calculation is
A_times_x = A_LT*x+(x.'*A_LT).'+ diag(A_LT).*x;
This should only perform transposes on the smaller resultant matrices.
7 个评论
Bruno Luong
2021-4-24
编辑:Bruno Luong
2021-4-24
Might be you can experiment with this as well
reshape((x.'*A),[],1)
Clayton Gotberg
2021-4-24
I begin to understand why engineers are specifically trained and hired to test software. It's deeply interesting to get this peek into what MATLAB does to ensure I can keep writing ill-considered code.
另请参阅
类别
在 Help Center 和 File Exchange 中查找有关 Logical 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!