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$
$\endgroup$
5
2 Answers
$\begingroup$
$\endgroup$
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.
answered Feb 5, 2020 at 7:26
$\begingroup$
$\endgroup$
4
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), becausemis the dimension of the Krylov subspace, andnis the system size, right? $\endgroup$– HappyFeb 5, 2020 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$ Feb 5, 2020 at 13:19
-
$\begingroup$ Thanks for the correction, it should have been n>>m, now updated. $\endgroup$ Feb 5, 2020 at 14:51
-
$\begingroup$ @FedericoPoloni, an interesting idea that I never heard.
mhas 3 legs, :) I like this. Thanks very much. $\endgroup$– HappyFeb 7, 2020 at 3:10
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$zeros(m,m), then it stores all the memory. And in fact, this is the more common way that we use, because most people usesmatlabto 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 calledzeros(m,m), because he usesmatlabto 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 callzeros(m,m), i.e., we cannot choosematlabto do this, right? $\endgroup$