This article presents a basic technique for singular value decomposition (SVD). It starts by presenting the One-Sided-Jacobi-Rotation for SVD of square (but not necessarily symmetric) matrices. It then shows that this method can be used to obtain the SVD of rectangular matrices, if a QR decomposition is carried out before the one sided Jacobi method is used to calculate the SVD of the R part of the QR decomposition. The article finishes by presenting some numerical outputs.
The following reference was useful for the algorithms
As usual, the code presented here is for study purposes only and should not be used in production systems.
The singular value decomposition (SVD) of a matrix \(A \in \mathbb{R}^{m \times n}\) is given by:
$$ A = U \Sigma V^{T} $$
where \(U \in \mathbb{R}^{m \times n}\) and \(V \in \mathbb{R}^{n \times n}\) are unitary (orthogonal) matrices, and \(\Sigma \in \mathbb{R}^{n \times n}\) is a diagonal matrix of singular values. If \(A\) is square symmetric, its SVD is the same as its eigenvalue decomposition. Throughout this article, we shall assume that the matrix \(A\) is of full rank.
The one-sided-Jacobi SVD is very similar to the Classic Jacobi eigenvalue decomposition given in a previous article except that the algorithm is applied to \(A^{T}A\) rather than \(A\). The inner part of the algorithm is given below:
\(proc \textit{ } One-Sided-Jacobi-Rotation (G, j, k)\)
\(\quad \textit{Compute } a_{jj} = \left(G^{T}G\right)_{jj}, a_{jk} = \left(G^{T}G\right)_{jk}, \textit{ and } a_{kk} = \left(G^{T}G\right)_{kk}\)
\(\quad \textit{if } |a_{jk}| \textit{is not too small} \)
\(\quad \quad \tau = (a_{jj} - a_{kk})/(2 \cdot a_{jk})\)
\(\quad \quad t = sign(\tau)/(|\tau| + \sqrt{1 + \tau^{2}})\)
\(\quad \quad c= 1/\sqrt{1 + t^{2}}\)
\(\quad \quad s = c \cdot t\)
\(\quad \quad G = G \cdot R(i, j, \theta) \quad \dots \textit{where } c = \cos{\theta} \textit{ and } s = \sin{\theta}\)
\(\quad \quad \textit{if right singular vectors are desired}\)
\(\quad \quad \quad J = J \cdot R(j, k, \theta) \)
\(\quad \quad \textit{end if} \)
\(\quad \textit{end if} \)
\(R(i, j, \theta)\) is the Givens rotation, see a previous article for more details. The SVD itself is given the in algorithm below:
\(\textit{repeat}\)
\(\quad \textit{for } j = 1 \text{ to } n - 1\)
\(\quad \quad \textit{for } k = j + 1 \text{ to } n\)
\(\quad \quad \quad \textit{call One-Sided-Jacobi-Rotation(G, j, k)}\)
\(\quad \quad \textit{end for}\)
\(\quad \textit{end for}\)
\(\textit{until } G^{T}G \textit{is diagonal enough}\)
\(\textit{Let } \sigma_{i} = \lVert G(:,i)\rVert_{2} (\textit{the 2-norm of column i of G})\)
\(\textit{Let } U = [u_{1}, \ldots, u_{n}], \textit{where } u_{i} = G(:,i)/\sigma_{i}\)
\(\textit{Let } V = J, \textit{the accumulated product of Jacobi rotations}\)
The QR SVD first carries out a QR decomposition on \(A = QR\) and then does a SVD on \(R = U_{1} \Sigma_{1} V_{1}\) using the above one-sided Jacobi method (or any other appropriate method), then the SVD of \(A\) is given as \(A = (G U_{1}) \Sigma_{1} V_{1}\).
The enumerations for statically dispatching the two methods is given below:
enum SVDMethod{OneSidedJacobi, JacobiQR}
alias OneSidedJacobi = SVDMethod.OneSidedJacobi;
alias JacobiQR = SVDMethod.JacobiQR;
The D code for implementing the one-sided Jacobi SVD is given below with informative comments. The implementation of the Jacobi rotation is similar to that given in a previous article
//Internal "workhorse" function for One-Sided Jacobi Rotation
private auto _one_sided_jacobi(T, T eps = max(1E-14, T.epsilon))(Matrix!(T) A)
{
immutable n = A.ncol;
//Left sided eigen vectors
auto U = zeros!(T)(A.nrow, n);
//Right sided eigen vectors
auto V = eye!(T)(n);
//For Eigen values
T[] sigma = new T[n];
T conv = T.infinity;
//Error threshold
immutable delta = eps*fnorm!(true)(gemm!(CblasTrans)(A, A));
do{
conv = 0;
for(long j = 0; j < n - 1; ++j)
{
for(long k = j + 1; k < n; ++k)
{
//Jacobi Rotation
auto R = SJacobi!(T, eps)(A, j, k);
//Apply the rotation
A = A.mult(R);
//Accumulate the rotation
V = V.mult(R);
conv += 2*pow(R.Apq, 2.0);
}
}
conv = conv.sqrt;
}while(conv > delta);
for(long i = 0; i < n; ++i)
{
sigma[i] = norm(A[i]);
for(long j = 0; j < n; ++j)
{
U[j, i] = A[j, i]/sigma[i];
}
}
auto obj = SVD!(T)(U.dup, sigma, V.dup);
return obj;
}
The SVD object is given by:
struct SVD(T)
if(isFloatingPoint!(T))
{
Matrix!(T) U;
T[] sigma;
Matrix!(T) V;
}
The One-Sided Jacobi QR is given by:
auto svd(SVDMethod method, T, T eps = max(1E-14, T.epsilon))(Matrix!(T) A)
if(method == OneSidedJacobi)
{
immutable n = A.ncol;
auto obj = _one_sided_jacobi!(T, eps)(A.dup);
auto idx = seq!(long)(0, n - 1);
sort!("a[0] > b[0]")(zip(obj.sigma, idx));
auto U = obj.U.dup;
auto V = obj.V.dup;
for(long i = 0; i < n; ++i)
{
obj.U[0..$, i] = U[0..$, idx[i]];
obj.V[0..$, i] = V[0..$, idx[i]];
}
return obj;
}
The implementation for the QR SVD is given below with comments:
auto svd(SVDMethod method, T, T eps = max(1E-14, T.epsilon))(Matrix!(T) A)
if(method == JacobiQR)
{
//QR Decomposition
auto qrA = qr!(ModifiedGramSchmidt)(A);
immutable n = A.ncol;
//Applying Jacobi Rotation to R
auto obj = _one_sided_jacobi!(T, eps)(qrA.R);
//Sorting according to eigen value size
auto idx = seq!(long)(0, n - 1);
sort!("a[0] > b[0]")(zip(obj.sigma, idx));
obj.U = gemm(qrA.Q, obj.U);
auto U = obj.U.dup;
auto V = obj.V.dup;
for(long i = 0; i < n; ++i)
{
obj.U[0..$, i] = U[0..$, idx[i]];
obj.V[0..$, i] = V[0..$, idx[i]];
}
return obj;
}
The SVD of a symmetric matrix A is given by:
auto A = matrix([2.326121, 1.711786, 2.502242, 1.715707, 1.417077, 2.319703, 1.851844,
1.711786, 2.219194, 2.163875, 1.932815, 1.467704, 2.551077, 1.899822,
2.502242, 2.163875, 3.473702, 2.084781, 2.135963, 3.192984, 2.467568,
1.715707, 1.932815, 2.084781, 1.998455, 1.626024, 2.465247, 1.820398,
1.417077, 1.467704, 2.135963, 1.626024, 1.753152, 2.243601, 1.686944,
2.319703, 2.551077, 3.192984, 2.465247, 2.243601, 3.717349, 2.735672,
1.851844, 1.899822, 2.467568, 1.820398, 1.686944, 2.735672, 2.392183], 7, 7);
"A:\n".writeln(A);
"SVD(A) (JacobiQR):\n".writeln(svd!(JacobiQR)(A));
output:
A:
Matrix(7 x 7)
2.326121 1.711786 2.502242 1.715707 1.417077 2.319703 1.851844
1.711786 2.219194 2.163875 1.932815 1.467704 2.551077 1.899822
2.502242 2.163875 3.473702 2.084781 2.135963 3.192984 2.467568
1.715707 1.932815 2.084781 1.998455 1.626024 2.465247 1.820398
1.417077 1.467704 2.135963 1.626024 1.753152 2.243601 1.686944
2.319703 2.551077 3.192984 2.465247 2.243601 3.717349 2.735672
1.851844 1.899822 2.467568 1.820398 1.686944 2.735672 2.392183
SVD(A) (JacobiQR):
SVD!double(Matrix(7 x 7)
0.3414 0.562 -0.5583 -0.06242 0.3344 0.2024 -0.3151
0.3433 -0.435 -0.4709 0.01698 -0.4495 -0.4122 -0.3146
0.4468 0.5458 0.2371 -0.07683 -0.486 -0.2401 0.3828
0.3348 -0.3483 -0.2316 -0.4309 0.3588 0.02034 0.631
0.3038 -0.06922 0.5325 -0.5201 0.2517 -0.1997 -0.4957
0.4761 -0.2563 0.1986 0.1966 -0.2241 0.7565 -0.08388
0.3673 -0.07289 0.1839 0.7036 0.4558 -0.3448 0.0649
, [15.4924, 0.895303, 0.632493, 0.40729, 0.265876, 0.161864, 0.0249189], Matrix(7 x 7)
0.3414 0.562 -0.5583 -0.06242 0.3344 0.2024 -0.3151
0.3433 -0.435 -0.4709 0.01698 -0.4495 -0.4122 -0.3146
0.4468 0.5458 0.2371 -0.07683 -0.486 -0.2401 0.3828
0.3348 -0.3483 -0.2316 -0.4309 0.3588 0.02034 0.631
0.3038 -0.06922 0.5325 -0.5201 0.2517 -0.1997 -0.4957
0.4761 -0.2563 0.1986 0.1966 -0.2241 0.7565 -0.08388
0.3673 -0.07289 0.1839 0.7036 0.4558 -0.3448 0.0649
)
We can compare this to the eigenvalue decomposition:
"eigen(A) (Pure):\n".writeln(eigen!(PureQR)(A));
output:
eigen(A) (Pure):
Eigen!double([15.4924, 0.895303, 0.632493, 0.40729, 0.265876, 0.161864, 0.0249189],
Matrix(7 x 7)
0.3414 -0.562 -0.5583 0.06242 0.3344 0.2024 -0.3151
0.3433 0.435 -0.4709 -0.01698 -0.4495 -0.4122 -0.3146
0.4468 -0.5458 0.2371 0.07683 -0.486 -0.2401 0.3828
0.3348 0.3483 -0.2316 0.4309 0.3588 0.02034 0.631
0.3038 0.06922 0.5325 0.5201 0.2517 -0.1997 -0.4957
0.4761 0.2563 0.1986 -0.1966 -0.2241 0.7565 -0.08388
0.3673 0.07289 0.1839 -0.7036 0.4558 -0.3448 0.0649
)
The SVD of an unsymmetric square matrix B is given by:
auto B = matrix([0.969015, 0.240537, 0.207134, 0.265424, 0.132492, 0.158843, 0.408621,
0.890505, 0.183104, 0.51322, 0.869608, 0.00461878, 0.704453, 0.640827,
0.431821, 0.264286, 0.660464, 0.859397, 0.973617, 0.78622, 0.161047,
0.935504, 0.237205, 0.849111, 0.623219, 0.276953, 0.21426, 0.636278,
0.881098, 0.0908641, 0.63171, 0.415055, 0.306431, 0.895314, 0.183507,
0.00495128, 0.197731, 0.934934, 0.940043, 0.731247, 0.31328, 0.104048,
0.717943, 0.787072, 0.362419, 0.320249, 0.930887, 0.68781, 0.528401], 7, 7);
"B:\n".writeln(B);
"SVD(B) (JacobiQR):\n".writeln(svd!(JacobiQR)(B));
output:
B:
Matrix(7 x 7)
0.969 0.8905 0.4318 0.9355 0.8811 0.004951 0.7179
0.2405 0.1831 0.2643 0.2372 0.09086 0.1977 0.7871
0.2071 0.5132 0.6605 0.8491 0.6317 0.9349 0.3624
0.2654 0.8696 0.8594 0.6232 0.4151 0.94 0.3202
0.1325 0.004619 0.9736 0.277 0.3064 0.7312 0.9309
0.1588 0.7045 0.7862 0.2143 0.8953 0.3133 0.6878
0.4086 0.6408 0.161 0.6363 0.1835 0.104 0.5284
SVD(B) (JacobiQR):
SVD!double(Matrix(7 x 7)
0.48 -0.6833 0.1013 0.01397 0.3628 0.383 0.118
0.206 0.0164 0.4826 -0.3458 -0.2667 -0.2969 0.6675
0.4248 0.2221 -0.4445 -0.1736 0.5154 -0.5217 0.06889
0.4405 0.2467 -0.5084 -0.03412 -0.5028 0.4381 0.2019
0.3515 0.5559 0.5071 -0.183 0.2207 0.3231 -0.3517
0.3955 0.04161 0.2082 0.8152 -0.1909 -0.3058 -0.06359
0.2698 -0.3345 0.01528 -0.3886 -0.44 -0.3212 -0.6061
, [3.7372, 1.32849, 0.953548, 0.691766, 0.428823, 0.258256, 0.0882725], Matrix(7 x 7)
0.2513 -0.454 0.09831 -0.2431 0.1862 0.6618 0.4351
0.4065 -0.3458 -0.3491 0.2237 -0.732 -0.07691 0.03852
0.4328 0.4427 0.1056 0.247 -0.02713 0.5471 -0.4945
0.3978 -0.2581 -0.3044 -0.5216 0.328 -0.2419 -0.4942
0.3757 -0.1593 -0.0148 0.6643 0.5324 -0.2749 0.1816
0.338 0.6204 -0.3774 -0.2624 0.03267 -0.09936 0.5278
0.413 0.03855 0.7887 -0.2182 -0.1914 -0.3359 0.09432
)
and the SVD of a rectangular matrix is given by:
auto C = matrix([0.5171, 0.720403, 0.886823, 0.545909, 0.191485, 0.164604,
0.525069, 0.0361945, 0.643514, 0.653749, 0.919446, 0.532787,
0.149449, 0.0414416, 0.846547, 0.0992408, 0.134474, 0.329745,
0.477403, 0.585913, 0.238699, 0.534982, 0.0995055, 0.140665,
0.0329689, 0.690778, 0.220595, 0.551619, 0.247787, 0.852184,
0.0890305, 0.692473, 0.397186, 0.526983, 0.793832, 0.386942,
0.791285, 0.253921, 0.981527, 0.949826, 0.81025, 0.448222,
0.656937, 0.014137, 0.197236, 0.45391, 0.615287, 0.646467,
0.84736, 0.611963], 10, 5);
"C:\n".writeln(C);
"SVD(C):\n".writeln(svd!(JacobiQR)(C));
output:
C:
Matrix(10 x 5)
0.5171 0.919446 0.238699 0.0890305 0.81025
0.720403 0.532787 0.534982 0.692473 0.448222
0.886823 0.149449 0.0995055 0.397186 0.656937
0.545909 0.0414416 0.140665 0.526983 0.014137
0.191485 0.846547 0.0329689 0.793832 0.197236
0.164604 0.0992408 0.690778 0.386942 0.45391
0.525069 0.134474 0.220595 0.791285 0.615287
0.0361945 0.329745 0.551619 0.253921 0.646467
0.643514 0.477403 0.247787 0.981527 0.84736
0.653749 0.585913 0.852184 0.949826 0.611963
SVD(C):
SVD!double(Matrix(10 x 5)
0.3136 0.6762 -0.2704 0.3668 -0.192
0.368 -0.06804 0.008703 -0.0641 -0.4595
0.2905 -0.2386 -0.1389 0.6305 -0.07492
0.1703 -0.4487 -0.1006 -0.04504 -0.3389
0.2689 0.1485 -0.5436 -0.6036 0.1124
0.2181 0.04894 0.5729 -0.09303 0.01993
0.3088 -0.3081 0.01678 0.08765 0.4326
0.2199 0.3724 0.4096 -0.01471 0.2813
0.425 -0.1422 -0.1542 0.04701 0.5106
0.4568 -0.03997 0.286 -0.2804 -0.3025
, [3.54787, 0.974648, 0.895078, 0.831674, 0.523949], Matrix(5 x 5)
0.4531 -0.3953 -0.2348 0.5285 -0.5513
0.3858 0.6688 -0.4927 -0.3122 -0.2521
0.3292 0.1488 0.8101 -0.2292 -0.4009
0.5459 -0.5213 -0.1101 -0.5469 0.3451
0.4893 0.3202 0.1837 0.5211 0.5939
)
Singular value decomposition has a very long history and the number of variants of the algorithm are numerous and highly involved. This article covered two very basic approaches. It is hoped that this topic will be revisited at a later stage to explore more algorithms for implementing this technique. Thank you for reading.