5
$\begingroup$

I have an $n \times n$ unsymmetric matrix $A$ that results from the discretization of an ill-posed Poisson problem, and thus is rank-deficient with null space of dimension one. I want to compute just the smallest singular triplet

\begin{equation*} (u_n,\sigma_n,v_n) \end{equation*}

where

\begin{equation*} u_n^T A v_n = \sigma_n = 0. \end{equation*}

I have found, through some experimentation, that applying inverse iteration to the shifted (and thus full-rank) matrix $A^T - \sigma I$ yields the left singular vector very quickly, often in one or two iterations. In Matlab code:

x(:,1) = rand(n,1); 
for k = 1:MAXIT
    x(:,k+1) = ( A' - sigma*eye(n) )\x(:,k); 
    x(:,k+1) = x(:,k+1)/norm(x(:,k+1)); 
end

Usually, this iteration converges to small enough tolerances ($|| x_k^T A|| = 10^{-9}$ or so) in one or two iterations.

I think this converges so rapidly because the largest and second largest singular values of $(A^T - \sigma I)^{-1}$ are spaced by $\sigma^{-1}$ which can be made large, but I have not been able to prove that this inverse iteration should converge. I'm sure there is some clever way of using the SVD to show this convergence as in the standard proof of the equivalent inverse iteration for the eigenvectors, but I can't come up with it.

Any help, ideas, or thoughts on why this works, and how I could prove this works would be helpful.

$\endgroup$
  • $\begingroup$ What do you tend to use for $\sigma$? Inverse iteration should converge to the eigenvector associated with the eigenvalue closest to $\sigma$. Apart from that, if you know $\sigma_n = 0$, then the left singular vector should be the eigenvector associated with the zero eigenvalue of $A^*$. Inverse iteration converges with rate $|\sigma - \lambda_{\text{closest to }\sigma}|/|\sigma-\lambda_{\text{next closest to }\sigma}|$. Can you say anything about the eigenvalues of $A^*$ (or equivalently, $A$)? $\endgroup$ – Jesse Chan Jan 6 '14 at 15:25
  • $\begingroup$ What method do you use to compute the inverse matrix-vector product in the inverse power iteration? If it is anything pivoted, then applying the same pivoted factorization method to the unmodified matrix $A$ should give you the kernel vectors using the last row and column of the factorization. $\endgroup$ – LutzL Jan 6 '14 at 17:27
  • 2
    $\begingroup$ $v_n$ is a kernel vector of $A$ of norm $1$, the same for $u_n$ and $A^\top$. If you compute $A=PLU$, then $Ax=0$ is equivalent to $Ux=0$ and by assumption and construction, the last row of $U$ is zero, which allows to set $x_{n}=1$ and solve the other equations in that triangular system. Then $v_n=x/\|x\|$. $\endgroup$ – LutzL Jan 6 '14 at 19:49
  • 1
    $\begingroup$ Pivoting in the QR decomposition requires column switches, so $AP=QR$ is the result, with $R$ again upper triangular and last diagonal entry equal to zero. Thus again, solve $Rx=0$ with $x_n=1$ and set $v_n=Px/\|x\|$. The left singular eigenvector is even easier to compute, it is $u_n=y/\|y\|$ with $y=Q^\top e_n$. $\endgroup$ – LutzL Jan 6 '14 at 19:56
  • 1
    $\begingroup$ Sumedh, I meant that if you have a zero singular value, this implies you have a zero eigenvalue as well - if $Ax = 0$, $A^*Ax = 0$ as well. Thus, your inverse iteration on $A^*$ is simply converging to the zero eigenvector, which is why it works. LutzL notes this as well in his PLU comments. $\endgroup$ – Jesse Chan Jan 6 '14 at 20:28
3
$\begingroup$

Just collecting the answers contained in the comments:

The left and right singular vectors of $A$ are exactly the left and right kernel vectors of $A$. Since the assumption is that the kernels have dimension $1$ and $A$ is well conditioned outside the kernel, one can compute the kernel vectors using the usual, pivoted, matrix factorizations.

For the LU factorization, $A=PLU$, $P$ a permutation matrix, $L$ lower triangular with $1$-diagonal, $Ax=0$ is equivalent to $Ux=0$. The last row of $U$ including the diagonal element is zero, thus choosing $x_n=1$ and substituting backwards results in a kernel vector $x$ and $v_n=x/\|x\|$.

$y^TA=0$ is equivalent to $(y^TPL)R=0$, which has the nontrivial solution $(y^TPL)=e_n^T$ or $L^T(P^Ty)=e_n$. $L^T$ is a regular triangular matrix, thus this is easily solvable and $y=P(L^{-T}e_n)$.


In the QR decomposition, row pivoting is without effect, so column pivoting needs to be applied, the factorization is then $AP=QR$. $Ax=0$ is equivalent to $R(P^Tx)=0$ which can be solved as before for $P^Tx$, and $y^TA=0$ is equivalent to $(y^TQ)R=0$, again with the obvious solution $y^TQ=e_n^T$ or $y=Qe_n$

$\endgroup$
  • $\begingroup$ Thanks for summarizing this for future readers. As a matter of courtesy, is it the role of the question asker to do this? I'm just getting my feet wet in scicomp, and don't know the rules of the road yet :) $\endgroup$ – Sumedh Joshi Jan 7 '14 at 16:11
  • $\begingroup$ No, it is the task of the responders to decide when a comment turns from requesting corrections or more information on the problem formulation to providing partial or full answers. That is, I should have put my second comment into an answer from the start, but it was not clear to me if there were other reasons to avoid direct factorization. $\endgroup$ – LutzL Jan 7 '14 at 16:43
  • $\begingroup$ I think you need complete pivoting instead of partial pivoting here. Otherwise the $LU$ factorization fails for $A = \begin{pmatrix} 0 & 2 \\ 0 & 1 \end{pmatrix}$. $\endgroup$ – wim Apr 20 at 8:19
  • $\begingroup$ .... Unless the $LU$ factorization code simply neglects the zero column, which seems to be the case in Octave. $\endgroup$ – wim Apr 22 at 8:15

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.