Hypothetical question for some code that I'm writing. Suppose I have an matrix-free linear operator $A$, i.e. the only thing I know about it is the forward action $v \mapsto Av$. For simplicity, let's further assume that $A$ is symmetric positive definite.

Are there any methods that would allow me to generate a 'decomposition' of this operator without needing to build it explicitly? Considering, say, a Cholesky decomposition: Is there a way to obtain an operator $L$ such that $L L^T v = Av$? Whether $L$ itself is a fully constructed matrix or not.

I suppose you could just use a regular Cholesky algorithm, and everytime you want to access element $A_{ij}$ you compute $e_i A e_j$, but perhaps there's a more performant or iterative way?

Stepping away from Cholesky, is there any 'factorization' I can realistically (meaning with just matrix-applies) perform on a matrix-free operator?

Thanks for any advice.

  • 2
    $\begingroup$ I can't comment specifically on computing a matrix-free computation of a Cholesky factor, but one can compute many matrix functions $f(A)$ using Krylov methods. $f(z) = z^{1/2}$ provides one such decomposition. Some resources can be found here scicomp.stackexchange.com/questions/11369/… $\endgroup$
    – whpowell96
    Commented May 23, 2023 at 23:00
  • 1
    $\begingroup$ @whpowell96 True, but it all boils down to what OP needs the Cholesky factor for. If it's to solve linear systems, then having $A^{1/2}$ instead of the Cholesky factor does not work; if it is for normalization purposes then the two might be equivalent. OP, if you could tell us more it would be helfpul; this might be an XY problem. $\endgroup$ Commented May 24, 2023 at 11:58
  • 1
    $\begingroup$ I think OP is missing the point of the matrix-free approach to begin with: rather than computing a Cholesky factor (or whatever else) and storing it in memory, we should ask ourselves what we would do with the Cholesky factor, and try to get at that directly using just matrix-vector products. Maybe I want a Cholesky factor because I want to solve a linear system as well as compute a determinant: then do those things directly using, say Lanczos. $\endgroup$ Commented May 24, 2023 at 19:23

1 Answer 1


This is not efficient, primarily because (i) the factors of sparse matrices are, in general, not sparse themselves, and (ii) when computing matrix factorizations, you don't just need $A_{ij}$ but you start with them and then you update them -- i.e., you have to allocate storage for them anyway, and so it isn't useful to try and start with a matrix-free representation: You need a matrix representation anyway.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.