2
$\begingroup$

Whenever I work with matrices for the SVD decomposition in the following manner:

X = np.random.randn(100,20)

U,S,V = np.linalg.svd(X)

plt.imshow(U) plt.figure() plt.imshow(V)

The larger matrix always exhibits this diagonally dominant behavior across the larger matrix (in this case $U$) as follows:

enter image description here

Currently I am not sure why this tends to appear at all. If the matrices were of a similar size such as if I did:

X = np.random.randn(100,95)

Then no diagonal component is observed in either $U$ or $V$. I am wondering why this is so.

Is this due to an artifact of the algorithm chosen by Python (in this case), or is there a strong mathematical reason to believe this always appears regardless of algorithmic method chosen to compute the SVD? Or perhaps is it a combination of both?

$\endgroup$
4
  • 1
    $\begingroup$ The trailing columns are not uniquely determined, and can be chosen to be any orthonormal basis that spans the same space. So it is an artifact of the algorithm used. Numpy's svd is a wrapper for lapack's dgesdd routine. On line 616 of dgesdd.f, the comment says that if the matrix is significantly taller than it is wide, it is first reduced in dimension via QR factorization. The QR routine used is dgeqrf, which does QR via Householder reflectors. So, the question reduces to asking why Householder QR on a tall matrix produce the diagonal structure. However I do not know the answer to that. $\endgroup$
    – Nick Alger
    Apr 28, 2023 at 4:39
  • $\begingroup$ Thanks for the comment. I came up with the same conclusion, that it is feels decomposition method driven (hence why I mentioned Python explicitly). But I don't know for example if this kind of structure appears quite universally regardless of the decomp method used (it certainly appeared in Python, for both taller and wider matrices). And if so, why? $\endgroup$ Apr 28, 2023 at 5:24
  • 1
    $\begingroup$ Here is some code showing similar results for Householder QR only: import numpy as np import matplotlib.pyplot as plt Q, R = np.linalg.qr(np.random.randn(100,20), mode='complete') plt.matshow(Q) $\endgroup$
    – Nick Alger
    Apr 28, 2023 at 21:19
  • 1
    $\begingroup$ After thinking a bit more about why Householder QR would create this phenomenon, I managed to convince myself that it makes sense. I wrote this up in an answer below. $\endgroup$
    – Nick Alger
    May 3, 2023 at 22:08

1 Answer 1

1
+50
$\begingroup$

The trailing columns may be any orthonormal basis that spans the orthogonal complement of the column space of $X$. Hence it is an artifact of the algorithm used. The question now becomes: what in the algorithm causes this diagonal behavior, and why does it stop occurring when the number of rows and columns are similar?

I did a little detective work about what algorithms numpy and matlab use, and in summary, the diagonal behavior is caused by a Householder QR decomposition that is performed as a pre-processing step if the matrix is significantly taller than it is wide. This is before the main SVD calculations begin.

Numpy and Matlab's svd code both wrap the dgesdd Fortran routine in LAPACK. On lines 616-618 of dgesdd.f, the comment says that if the matrix is significantly taller than it is wide, it is first reduced in dimension via full QR factorization. Here I quote lines 614-635:

dgesdd.f
...
614       IF( m.GE.n ) THEN
615 *
616 *        A has at least as many rows as columns. If A has sufficiently
617 *        more rows than columns, first reduce using the QR
618 *        decomposition (if sufficient workspace available)
619 *
620          IF( m.GE.mnthr ) THEN
621 *
622             IF( wntqn ) THEN
623 *
624 *              Path 1 (M >> N, JOBZ='N')
625 *              No singular vectors to be computed
626 *
627                itau = 1
628                nwork = itau + n
629 *
630 *              Compute A=Q*R
631 *              Workspace: need   N [tau] + N    [work]
632 *              Workspace: prefer N [tau] + N*NB [work]
633 *
634                CALL dgeqrf( m, n, a, lda, work( itau ), work( nwork ),
635      $                      lwork - nwork + 1, ierr )
...

Mathematically, what they are doing is first factoring $$X = QR = \begin{bmatrix}Q_1 & Q_2\end{bmatrix}\begin{bmatrix}R_1 \\ 0\end{bmatrix}$$ where $Q$ is orthonormal, and $R_1$ is square and upper triangular. Then the SVD only needs to be performed on the smaller matrix $R_1$: $$R_1 = U' \Sigma' V'^T$$ from which they recover the SVD of $X$ as follows: $$X = U \Sigma V^T = \begin{bmatrix}Q_1 U' & Q_2\end{bmatrix} \begin{bmatrix}\Sigma' \\ 0\end{bmatrix} V'^T.$$ The diagonal behavior you are seeing is in the matrix $Q_2$.

The QR routine used is dgeqrf, which does QR via Householder reflectors. This works as follows. First, define $$u := x - \text{sign}(x_0)||x||e_1, \qquad v := u / ||u||, \qquad Q^{(1)} := I - 2 v v^T$$ where $x$ is the first column of $X$, $\text{sign}(x_0)$ is the sign of the first entry of $x$, and $e_1=(1,0,0,\dots,0)$. Then the action of $Q^{(1)}$ zeros out the first column of $X$: $$Q^{(1)} X = \begin{bmatrix}||x|| & \text{stuff} \\ 0 & X'\end{bmatrix}.$$ Because of how $Q^{(1)}$ is constructed from $x$, the first column and row of $Q^{(1)}$ will be a normalized version of $x$, while the rest of the matrix $Q$ will be the identity plus the outer product of a random vector with itself. In other words, a perturbation of the identity by small noise. Here is some code and a picture for a small example:

N = 10
x = np.random.randn(N)
e0 = np.zeros(N)
e0[0] = 1.0
a = -np.sign(x[0]) * np.linalg.norm(x)
u = x - a * e0
v = u / np.linalg.norm(u)
Q = np.eye(N) - 2.0 * np.outer(v, v)
plt.matshow(Q)

QR householder first Q

The process then repeats on the smaller matrix $X'$, and so on, until $X$ has been transformed into $$R = Q^{(k)} Q^{(k-1)} \dots Q^{(1)} X$$ by actions of a sequence of orthogonal matrices of the form $$Q^{(i)} = \begin{bmatrix}I & 0 \\ 0 & I - 2 v v^T\end{bmatrix}.$$ Here the $v$ is the analogous vector for the smaller submatrix that is being worked on in the ith step. By multiplying these matrices together to form $Q = Q^{(k)} Q^{(k-1)} \dots Q^{(1)}$, we end up simply adding random noise repeatedly to the identity in the bottom right submatrix, which yields the diagonal behavior observed.

$\endgroup$

You must log in to answer this question.

Not the answer you're looking for? Browse other questions tagged .