Suppose one has a large sparse symmetric positive definite matrix $A$ and wants to multiply it by a vector $x$. Only the lower triangular part of matrix A is stored/known. The multiplication $Ax$ should be as efficient as possible, as it has to be performed many times and $A$ is large (millions of rows and columns). Is there any known algorithm that achieves this?
-
$\begingroup$ math.stackexchange.com/q/4118934/14578 $\endgroup$– D.W. ♦Jun 8, 2021 at 5:36
-
$\begingroup$ It is my understanding that Matlab-specific questions are off-topic here. $\endgroup$– D.W. ♦Jun 8, 2021 at 5:37
-
$\begingroup$ @D.W. ok, this is a question regarding algorithms. If it happens that anyone knows about a matlab implementation I appreciate it. Thank you $\endgroup$– MarxJun 8, 2021 at 11:04
-
$\begingroup$ How is the sparse matrix stored precisely ? $\endgroup$– Yves DaoustMar 17, 2022 at 9:25
2 Answers
Answer by Clayton Gotberg [1], modified:
If $\textbf{A}$ is a symmetric matrix and $\textbf{A}_{LT}$ is the lower triangular part of the matrix and $\textbf{A}_{UT}$ is the upper triangular part of the matrix:
$\textbf{A}_x = \textbf{A}_{LT}x + \textbf{A}_{UT}x - diag(\textbf{A}) \cdot x$
where the diagonal function only finds the diagonal elements of $\textbf{A}$. This is because of a few relations:
$\textbf{A} = \textbf{A}_{LT} + \textbf{A}_{UT} - \textbf{IA}$ $(\textbf{B} + \textbf{C})x = \textbf{B}x + \textbf{C}x$
To save time and space [...], take advantage of the relations:
$\textbf{A}_{UT} = \textbf{A}_{LT}$ $(\textbf{A}x)^T = x^T\textbf{A}^T$
To get:
$\textbf{A}x = \textbf{A}_{LT}x + (x^T\textbf{A}_{LT})^T - diag(\textbf{A}_{LT}) \cdot x$
Additionally, there are efficient GPU implementations [2] for efficient triangular matrix vector multiplication - this paper includes several references for relevant papers and algorithms.
References
[1]: Symmetric Matrix-Vector Multiplication with only lower triangular stored, https://au.mathworks.com/matlabcentral/answers/812175-symmetric-matrix-vector-multiplication-with-only-lower-triangular-stored/?s_tid=ans_lp_feed_leaf
[2]: Efficient Triangular Matrix Vector Multiplication on the GPU, https://link.springer.com/chapter/10.1007%2F978-3-030-43229-4_42
As far as I know, there is no better way than the standard matrix-vector accumulation, taking $O(n^2)$ operations. Anyway, you need to take into account the packed storage of the matrix, and you need to process the matrix elements following an L shape.