For the Krylov subspace method to solve the large sparse linear system, we first need to generate a subspace Km = span{v,Av,...A^{m-1}v}, which indeed a process that we convert a linearly independent basis to an orthonormal basis: $$ V_mAV_m = H_m, $$ where $H_m$ is a upper Hessenberg matrix, which has about $m^2/2$ nonzero entries. So, when we code in matlab, we first need to H = zeros(m,m), so in fact, we need to store all the $m^2$ entries including the zero entries. But in many books, the author said the storage of the matrix $H_m$ is about $m^2/2$. I do not understand why author said that? Because the code H = zeros(m,m) indeed store $m^2$ entries even though $m^2/2$ zero entries. any suggestions?

  • $\begingroup$ You can of course use a data structure that really only stores the upper half of the matrix, at a memory cost of roughly $m^2/2$. The fact that you call zeros(m,m) and get a matrix that also has memory allocated for everything below the diagonal is your choice. It is not necessary for the algorithm. $\endgroup$ – Wolfgang Bangerth Feb 5 '20 at 2:56
  • $\begingroup$ @WolfgangBangerth, Thanks Prof. I get your comment. I know it is not necessary. But If one uses zeros(m,m), then it stores all the memory. And in fact, this is the more common way that we use, because most people uses matlab to code the algorithm. And I also found that some authors states that we need to store $m^2/2$ memory in his book, but he also called zeros(m,m), because he uses matlab to do the algorithm. So, I am sometimes a little confused. So, if we state that the memory is $m^2/2$, then we should not call zeros(m,m), i.e., we cannot choose matlab to do this, right? $\endgroup$ – sunshine Feb 5 '20 at 3:19
  • 3
    $\begingroup$ So then don't use matlab, or if you use matlab, use a better suited data structure. It's like saying "Some authors state that it takes 3 hours to get from Austin, TX to Houston TX, but when I ride my bike it takes all day", to which the answer is "Then don't take the bike, take a car!" Just because you don't use the right tool for a task doesn't invalidate the statement that it only takes $m^2/2$ memory :-) $\endgroup$ – Wolfgang Bangerth Feb 5 '20 at 13:08
  • 1
    $\begingroup$ @WolfgangBangerth, an interesting comment, Get it. it is just because we donot use the right tool, haha, thanks very much. ;) $\endgroup$ – sunshine Feb 7 '20 at 3:13
  • $\begingroup$ No, to be precise, it's because you don't use the right tool ;-) In professional implementations, we actually do often only store the nonzero entries ;-)) $\endgroup$ – Wolfgang Bangerth Feb 7 '20 at 4:15

Your view is correct, Matlab does store the zeros. As pointed out by @rchilton1980, this particular non-optimization that you are pointing out here is not too harmful, since the bulk of the storage in Krylov methods is the matrix V, not H. But that is just an instance of a general phenomenon. That is just Matlab's choice of a tradeoff between simplicity and efficiency.

Matlab often does not care about these sorts of optimizations. Another classical example is QR factorization, where it returns the full Q matrix instead of a compact representation that could use about $n^2/2$ storage. Or LU factorizations, where two separate matrices L and U are returned for simplicity. Factorizations are computed in new matrices and not in-place overwriting the inputs.

If you care about that level of efficiency, you will need to use Lapack routines directly, or switch to Julia which is better in this regard.


In practice I don't think that particular factor of $2$ is worth pursuing. Thinking of gmres applied to a large problem ($n >> m$), the $m^2$ storage of that small Hessenberg projection is dwarfed by the $m\cdot n$ storage of the search vectors themselves (the $\mathbf V_m$). You are forced into restarts because $\mathbf V_m$ is too big, not because of $\mathbf H_m$.

  • 1
    $\begingroup$ Thanks very much. a question is that in your comments, a large problems (m>>n), it seems wrong, it should be (m < < n), because m is the dimension of the Krylov subspace, and n is the system size, right? $\endgroup$ – sunshine Feb 5 '20 at 2:28
  • 1
    $\begingroup$ Usually, I like the convention that $m>n$ because $m$ has three legs and $n$ has two. And because it is consistent with the notation from graph theory that $m=|E|$, $n = |V|$. $\endgroup$ – Federico Poloni Feb 5 '20 at 13:19
  • $\begingroup$ Thanks for the correction, it should have been n>>m, now updated. $\endgroup$ – rchilton1980 Feb 5 '20 at 14:51
  • $\begingroup$ @FedericoPoloni, an interesting idea that I never heard. m has 3 legs, :) I like this. Thanks very much. $\endgroup$ – sunshine Feb 7 '20 at 3:10

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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