Singular Value Decomposition

For any matrix A ∈ ℝ^m×n there exist orthogonal matrices U ∈ ℝ^m×m, V ∈ ℝ^n×n and a ’diagonal’ matrix Σ ∈ ℝ^m×n with diagonal entries
σ₁ ≥ · · · ≥ σ_r > σ_r+1 = · · · = σ_min{m,n} = 0
such that A = UΣV^T.

The decomposition A = UΣV^T is called Singular Value Decomposition(SVD). It has fundamental importance decomposition of a matrix and tells us a lot about its structure. It is very important in several linear algebra applications.
The diagonal entries σ_i of Σ are called the singular values of A. The columns of U are called left singular vectors and the columns of V are called right singular vectors.
Using the orthogonality of V we can write it in the form AV = UΣ. We can interpret it as follows: there exists a special orthonormal set of vectors (i.e. the columns of V ), that is mapped by the matrix A into an orthonormal set of vectors (i.e. the columns of U).

In core.matrix, the linear/svd function performs SVD and returns a map containing the key-value pairs [:U :V* :S], where :U and :V(V* is same as the transpose of V, i.e. V^T) are the orthogonal matrices containing the left and right singular vectors and :S is a vector of all the singular values of A(generally in decreasing order).
If you are only interested in a particular matrix and would like to avoid wastage of space and time, the API allows you to specify which matrices to return by passing a :return parameter in the options map.

Intended usage:

(let [{:keys [U S V*]} (svd M)] ....)
; => returns a map containing all three results, U, V* and S
(let [{:keys [S]} (svd M {:return [:S]})] ....)
; => returns a map containing just S

Given an SVD decomposition A = UΣV^T, σ₁ ≥ · · · ≥ σ_r > σ_r+1 = · · · = σ_min{m,n} = 0, if we denote columns of U as u₁, ..., u_m and columns of V as v₁, . . ., v_n, we have:

• rank(A) = r
• R(A) = R([u1, . . . , ur])
• N(A) = R([vr+1, . . . , vn])
• R(A^T) = R([v1, . . . , vr])
• N(A^T) = R([ur+1, . . . , um])

Also, if U_r = [u₁, . . . , u_r], Σ_r = diag(σ₁, . . . , σ_r), V_r = [v₁, . . . , v_r], then we have,
A = U_rΣ_rV^T_r = Σ^r_i=1(σ_iu_iv_i^T).
This is called the dyadic decomposition of A, decomposes the matrix A of rank r into sum of r matrices of rank 1.

Frobenius norm of a matrix can be easily computed from the SVD decomposition:
||A||_F = Σ^m_i=1Σⁿ_j=1a²_ij = sqrt(σ₁² + · · · + σ_p²), p = min {m, n}

From the SVD decomposition of A, it also follows that A^TA = VΣ^TΣV^T and AA^T = UΣ^TΣU^T
Thus, σ²_i, i = 1, . . . , p are the eigenvalues of symmetric matrices A^TA and AA^T and v_i and u_i are the corresponding eigenvectors.

Another important application of SVD is solving the linear least squares problem. Given A^{m, n} and b ∈ ℝ^m, with m ≥ n ≥ 1. The problem to find x ∈ ℝⁿ that minimizes ||Ax-b||₂ is called the least squares problem.
A minimizing vector x is called the least squares solution of Ax = b.
let A = UΣV^T be the SVD of A ∈ ℝ^m×n
||Ax-b||²₂ = ||AVV^Tx-b||²₂                        (VV^T = I)
                = ||U^T(AVV^Tx-b)||²₂                 (Orthogonal matrices are isometric)
                = ||ΣV^Tx-U^Tb||²₂                      (A = UΣV^T)
let V^Tx be z
                = Σ^r_i=0(σⁱzⁱ - u^Tib)² + Σ^m_i=r+1(u^T_ib)²

Thus we need to find x that for which the following is minimum
                Σ^r_i=0(σⁱzⁱ - u^Tib)² + Σ^m_i=r+1(u^T_ib)²
Therefore,
                z_i = u_i^Tb/σ_i, i=1,...r and
                z_i = artitrary, i=r+1,...n
As a result, x is minimum when
                ||Ax-b||²₂ = Σ^m_i=r+1(u^T_ib)²

Since z = V^Tx and V is orthogonal
                ||x||₂ = ||VV^Tx||₂ = ||V^Tx||2 = ||z||₂

The minimum norm solution of the linear least squares problem is given by
                z = V^Tx,
i.e.            x = Vz,

Therefore, the minimum norm solution is
                x = Σ^r_i=0(u_i^Tbv_i/σⁱ)