Before I start: The technique presented in this post is only an investigation to compute SVD. The algorithm presented here is (possibly) not numerically stable, and is (probably) not suitable when high accuracy (beyond 8 decimals) is needed. And the algorithm is very slow compared to standard SVD algorithms. That said, for small matrices, the technique for SVD presented here is the simplest I’ve seen (but by no means simple) and appears to work well.
The problem is singular value decomposition (SVD). If you have a matrix A and decompose it using SVD, you get a matrix U, a vector s, and a matrix Vh, where A = U * S * Vh. The S is a square matrix with the s values on the diagonal. SVD has many important uses in machine learning, for example, computing a pseudo-inverse for training a linear model.
Unfortunately, accurately computing SVD is incredibly difficult. Many mathematicians worked for many years on SVD and there are dozens of solutions, each with pros and cons. Two common SVD algorithms are Jacobi and Golub-Kahan-Reinsch. Both algorithms are very complex. But it is possible to compute an approximation to SVD. Briefly, to approximate the SVD of a source matrix A, first compute the QR decomposition of A to get Q and R, and then compute the SVD of the R matrix.
This seems silly at first — using SVD to compute SVD — but the R result matrix of a QR decomposition is a square matrix, and it’s possible to find its SVD using QR in a relatively simple way . In high-level pseudo-code:
def my_svd(A): Q, R = my_qr(A, mode='reduced') # apply QR decomp U, s, Vh = my_svd_square(R) # apply QR decomp again U = np.matmul(Q, U) return U, s, Vh
The SVD function described above isn’t trivial, but it’s vastly less complicated than computing SVD using Jacobi or any other algorithm. The my_svd_square() helper function for a square matrix is also much simpler than computing SVD for a non-square matrix. In essence, the idea for approximating SVD is to break it down into two much simpler (but by no means easy) sub-problems that can be solved using QR decomposition.
A few days ago, I implemented the idea using Python. Just for fun, one rainy morning before work, I decided to refactor the demo code to the C# language. The output of a demo run:
Begin SVD using QR demo with C# Source matrix A: 1.0000 6.0000 8.0000 5.0000 0.0000 2.0000 4.0000 6.0000 9.0000 3.0000 7.0000 7.0000 -1.0000 -3.0000 -5.0000 -7.0000 7.0000 5.0000 3.0000 1.0000 Computing SVD using QR Exited at iteration 20 U = 0.490569 0.357946 -0.658897 -0.443928 0.309393 0.418534 0.228973 0.339519 0.618658 -0.402788 0.526177 -0.422092 -0.402797 -0.389987 -0.198233 -0.465342 0.344434 -0.618365 -0.444148 0.541249 s = 21.095067 8.254694 4.718203 1.611445 Vh = 0.420588 0.395767 0.594453 0.559554 -0.872923 -0.017625 0.219636 0.435263 0.247105 -0.750910 -0.214769 0.573539 0.007029 0.528386 -0.743142 0.410485 Checking U * S * Vh: 1.0000 6.0000 8.0000 5.0000 0.0000 2.0000 4.0000 6.0000 9.0000 3.0000 7.0000 7.0000 -1.0000 -3.0000 -5.0000 -7.0000 7.0000 5.0000 3.0000 1.0000 End demo
I verified my results using the NumPy library np.linalg.svd() function and got the same results — essentially. One of the hugely annoying details with SVD is that the sign of the values in the columns of U, and in the rows of Vh are arbitrary. Singular values by convention are always positive or zero.
Well, all in all, this was a very challenging and very interesting exploration.
The history of machine learning is filled with incremental innovations. The history of pinball machines is filled with incremental innovations. A major advance occurred in 1953 when “Army Navy” by Williams introduced the first automatic reel scoring. Before this, fixed scores (1000 2000 etc.) were printed on the back glass and were illuminated as a player reached a point threshold. Somewhat surprisingly, the reel scoring feature did not catch on until a few years later when a high score reel was added.
Demo program. Replace “lt” (less than), “gt”, “lte”, “gte” with Boolean operator symbols — my blog editor chokes on symbols.
using System;
using System.IO;
using System.Collections.Generic;
namespace Matrix_SVD_Using_QR
{
internal class Program
{
static void Main(string[] args)
{
Console.WriteLine("\nBegin SVD using QR" +
" demo with C# ");
double[][] A = new double[5][];
A[0] = new double[] { 1, 6, 8, 5 };
A[1] = new double[] { 0, 2, 4, 6 };
A[2] = new double[] { 9, 3, 7, 7 };
A[3] = new double[] { -1, -3, -5, -7 };
A[4] = new double[] { 7, 5, 3, 1 };
Console.WriteLine("\nSource matrix A: ");
MatShow(A, 4, 9);
Console.WriteLine("\nComputing SVD using QR ");
double[][] U; double[] s; double[][] Vh;
MatSVD(A, out U, out s, out Vh);
Console.WriteLine("\nU = ");
MatShow(U, 6, 11);
Console.WriteLine("\ns = ");
VecShow(s, 6, 11);
Console.WriteLine("\nVh = ");
MatShow(Vh, 6, 11);
Console.WriteLine("\nChecking U * S * Vh: ");
double[][] S = MatMake(s.Length, s.Length);
for (int i = 0; i "lt" s.Length; ++i)
S[i][i] = s[i];
double[][] C = MatProduct(U, MatProduct(S, Vh));
MatShow(C, 4, 9);
Console.WriteLine("\nEnd demo ");
Console.ReadLine();
}
// ------------------------------------------------------
private static void MatSVD(double[][] A,
out double[][] U, out double[] s, out double[][] Vh)
{
// an approximation for SVD using QR
// not practical for anything other than small matrices
double[][] Q; double[][] R;
MatDecomposeQR(A, out Q, out R, reduced:true);
double[][] UU; double[] ss; double[][] VVh; // working
MatSVD_Square(R, out UU, out ss, out VVh);
UU = MatProduct(Q, UU);
U = UU; s = ss; Vh = VVh;
}
// ------------------------------------------------------
private static void MatSVD_Square(double[][] A,
out double[][] U, out double[] s, out double[][] Vh)
{
// SVD for a square matrix using iterative QR
int m = A.Length; int n = A[0].Length;
if (m != n)
throw new Exception("MatSVD_Square accepts" +
" only square matrix ");
int maxIter = 1000;
double[][] R1 = MatCopy(A);
double[][] UU = MatIdentity(m);
double[][] VV = MatIdentity(m);
double[][] Q1; double[][] Q2; double[][] R2;
for (int i = 0; i "lt" maxIter; ++i)
{
MatDecomposeQR(R1, out Q1, out R1, reduced:true);
MatDecomposeQR(MatTranspose(R1),
out Q2, out R2, reduced:true);
R1 = MatTranspose(R2);
UU = MatProduct(UU, Q1);
VV = MatProduct(VV, Q2);
// exit when off-diag cells of R1 are all close to 0
if (i % 10 == 0)
{
if (MatOffDiagAllZero(R1, 1.0e-8) == true)
{
Console.WriteLine("Exited at iteration " + i);
break;
}
}
} // i
// deal with negative singular values (in R1)
double[] ss = new double[m];
for (int j = 0; j "lt" m; ++j)
ss[j] = R1[j][j];
double[][] VVh = MatTranspose(VV);
for (int i = 0; i "lt" m; ++i)
{
if (ss[i] "lt" 0.0)
{
ss[i] = -ss[i]; // sing vals must be non-neg
// adjust row of VVh
for (int j = 0; j "lt" m; ++j)
VVh[i][j] = -VVh[i][j];
}
}
U = UU; s = ss; Vh = VVh;
} // MatSVD_Square()
// ------------------------------------------------------
private static bool MatOffDiagAllZero(double[][] A,
double tol)
{
int nr = A.Length; int nc = A[0].Length;
for (int i = 0; i "lt" nr; ++i)
{
for (int j = 0; j "lt" nc; ++j)
{
if (i == j) continue; // no check diag
if (Math.Abs(A[i][j]) "gt" tol)
return false;
}
}
return true;
}
// ------------------------------------------------------
private static void MatDecomposeQR(double[][] M,
out double[][] Q, out double[][] R, bool reduced)
{
// QR decomposition, Householder algorithm.
// see rosettacode.org/wiki/QR_decomposition
int m = M.Length;
int n = M[0].Length;
if (m "lt" n)
throw new Exception("No rows less than cols");
double[][] QQ = MatIdentity(m); // working Q
double[][] RR = MatCopy(M); // working R
int end;
if (m == n) end = n - 1;
else end = n;
for (int i = 0; i "lt" end; ++i)
{
double[][] H = MatIdentity(m);
double[] a = new double[m - i]; // corr
int k = 0;
for (int ii = i; ii "lt" m; ++ii) // corr
a[k++] = RR[ii][i];
double normA = VecNorm(a);
if (a[0] "lt" 0.0 && normA "gt" 0.0) // corr
normA = -normA;
else if (a[0] "gt" 0.0 && normA "lt" 0.0)
normA = -normA;
double[] v = new double[a.Length];
for (int j = 0; j "lt" v.Length; ++j)
v[j] = a[j] / (a[0] + normA);
v[0] = 1.0;
// Householder algorithm
double[][] h = MatIdentity(a.Length);
double vvDot = VecDot(v, v);
double[][] A = VecToMat(v, v.Length, 1);
double[][] B = VecToMat(v, 1, v.Length);
double[][] AB = MatProduct(A, B);
for (int ii = 0; ii "lt" h.Length; ++ii)
for (int jj = 0; jj "lt" h[0].Length; ++jj)
h[ii][jj] -= (2.0 / vvDot) * AB[ii][jj];
// copy h[][] into lower right corner of H[][]
int d = m - h.Length; // corr
for (int ii = 0; ii "lt" h.Length; ++ii)
for (int jj = 0; jj "lt" h[0].Length; ++jj)
H[ii + d][jj + d] = h[ii][jj];
QQ = MatProduct(QQ, H);
RR = MatProduct(H, RR);
} // i
if (reduced == false)
{
Q = QQ; // working results into the out params
R = RR;
return;
}
//else if (reduced == true)
{
int qRows = QQ.Length; int qCols = QQ[0].Length;
int rRows = RR.Length; int rCols = RR[0].Length;
// assumes m "gte" n !!
// square-up R
int dim = Math.Min(rRows, rCols);
double[][] Rsquared = MatMake(dim, dim);
for (int i = 0; i "lt" dim; ++i)
for (int j = 0; j "lt" dim; ++j)
Rsquared[i][j] = RR[i][j];
// Q needs same number columns as R
// so that inv(R) * trans(Q) works
double[][] Qtrimmed = MatMake(qRows, dim);
for (int i = 0; i "lt" qRows; ++i)
for (int j = 0; j "lt" dim; ++j)
Qtrimmed[i][j] = QQ[i][j];
Q = Qtrimmed;
R = Rsquared;
return;
}
} // MatDecomposeQR()
// ------------------------------------------------------
private static double[][] MatCopy(double[][] M)
{
int nr = M.Length; int nc = M[0].Length;
double[][] result = MatMake(nr, nc);
for (int i = 0; i "lt" nr; ++i)
for (int j = 0; j "lt" nc; ++j)
result[i][j] = M[i][j];
return result;
}
// ------------------------------------------------------
private static double[][] MatTranspose(double[][] m)
{
int nr = m.Length;
int nc = m[0].Length;
double[][] result = MatMake(nc, nr); // note
for (int i = 0; i "lt" nr; ++i)
for (int j = 0; j "lt" nc; ++j)
result[j][i] = m[i][j];
return result;
}
// ------------------------------------------------------
private static double[][] MatMake(int nRows, int nCols)
{
double[][] result = new double[nRows][];
for (int i = 0; i "lt" nRows; ++i)
result[i] = new double[nCols];
return result;
}
// ------------------------------------------------------
private static double[][] MatIdentity(int n)
{
double[][] result = MatMake(n, n);
for (int i = 0; i "lt" n; ++i)
result[i][i] = 1.0;
return result;
}
// ------------------------------------------------------
private static double[][] MatProduct(double[][] matA,
double[][] matB)
{
int aRows = matA.Length;
int aCols = matA[0].Length;
int bRows = matB.Length;
int bCols = matB[0].Length;
if (aCols != bRows)
throw new Exception("Non-conformable matrices");
double[][] result = MatMake(aRows, bCols);
for (int i = 0; i "lt" aRows; ++i) // each row of A
for (int j = 0; j "lt" bCols; ++j) // each col of B
for (int k = 0; k "lt" aCols; ++k)
result[i][j] += matA[i][k] * matB[k][j];
return result;
}
// ------------------------------------------------------
private static double VecNorm(double[] vec)
{
int n = vec.Length;
double sum = 0.0;
for (int i = 0; i "lt" n; ++i)
sum += vec[i] * vec[i];
return Math.Sqrt(sum);
}
// ------------------------------------------------------
private static double[][] VecToMat(double[] vec,
int nRows, int nCols)
{
double[][] result = MatMake(nRows, nCols);
int k = 0;
for (int i = 0; i "lt" nRows; ++i)
for (int j = 0; j "lt" nCols; ++j)
result[i][j] = vec[k++];
return result;
}
// ------------------------------------------------------
private static double VecDot(double[] v1, double[] v2)
{
double result = 0.0;
int n = v1.Length;
for (int i = 0; i "lt" n; ++i)
result += v1[i] * v2[i];
return result;
}
// ------------------------------------------------------
private static void MatShow(double[][] M, int dec,
int wid)
{
for (int i = 0; i "lt" M.Length; ++i)
{
for (int j = 0; j "lt" M[0].Length; ++j)
{
double v = M[i][j];
Console.Write(v.ToString("F" + dec).
PadLeft(wid));
}
Console.WriteLine("");
}
}
// ------------------------------------------------------
private static void VecShow(double[] vec,
int dec, int wid)
{
for (int i = 0; i "lt" vec.Length; ++i)
Console.Write(vec[i].ToString("F" + dec).
PadLeft(wid));
Console.WriteLine("");
}
} // Program
} // ns
.NET Test Automation Recipes
Software Testing
SciPy Programming Succinctly
Keras Succinctly
R Programming
Visual Studio Live
Microsoft MLADS Conference
DevIntersection Conference
Machine Learning Week
Ai4 Conference
G2E Conference
iSC West Conference
You must be logged in to post a comment.