I enjoy writing code. I write code almost every day. One downside to writing lots of code is that I often refactor a technique several times using different algorithms and different design patterns — and so I end up with many different versions of a function or system.
For example, there are many ways to implement relaxed Moore-Penrose pseudo-inverse for machine learning scenarios to train a linear model (linear regression and quadratic regression are two common examples). For MP pseudo-inverse, there are several algorithms that use singular value decomposition (SVD), and several algorithms that use QR decomposition.
Based on my experience, implementing any of the SVD versions from scratch is not a good idea — SVD is one of the most complex algorithms in all of numerical programming. The standard LAPACK library implementation of SVD is over 10,000 lines of code. Implementing a QR decomposition version of MP pseudo-inverse from scratch is more feasible. There are three main QR decomposition algorithms: Householder, modified Gram-Schmidt, and Givens. Householder is the most complicated to implement but most robust in theory. Modified Gram-Schmidt is the easiest to implement but less robust than Householder in theory. The Givens algorithm is sort of in between — moderately complex to implement, and sort of moderately robust (but the robustness claim is somewhat dubious in my opinion).
One evening when I couldn’t sleep (which is every evening because I have a no-sleep gene inherited from my father), I decided to wade through my dozens of implementations of MP pseudo-inverse using QR decomposition and consolidate down to three version — one Householder, one Gram-Schmidt, one Givens.
I organized the code into three container classes: QRHouseholder, QRGramSchmidt, QRGivens. Each class exposes a MatPseudoInv() method which is computed using the associated QR decomposition algorithm. Here’s the Householder version but all three versions are identical; they just call different versions of MatDecompQR():
public static double[][] MatPseudoInv(double[][] M)
{
// relaxed Moore-Penrose pseudo-inverse
// using QR Householder decomposition
// A = Q*R, pinv(A) = inv(R) * trans(Q)
int nr = M.Length; int nc = M[0].Length; // aka m, n
if (nr < nc)
Console.WriteLine("ERROR: Works only m >= n");
double[][] Q; double[][] R;
MatDecompQR(M, out Q, out R); // Householder
double[][] Rinv = MatInvUpperTri(R); // std algo
double[][] Qinv = MatTranspose(Q); // is inv(Q)
double[][] result = MatProduct(Rinv, Qinv);
return result;
}
I tested all three versions of the relaxed Moore-Penrose pseudo-inverse functions via the three different QR decomposition algorithms. For each algorithm, I generated 5,000 different random matrices, with between 100 and 1,000 rows, and between 2 and 20 columns, and where each value is between -10.0 and +10.0. All three versions passed all tests. This isn’t a great set of tests by any means, but it suggests the three pseudo-inverse implementations are probably OK for non-critical work.
Most of the guys I work with enjoy all kinds of puzzles. I like the “Where’s Waldo” books, and I’ll bet that if you’re reading this blog post, you enjoy Waldo too.
The covers of Playboy Magazine are sort of like Where’s Waldo for adults. Every cover, except for the very first one in 1953, has a bunny logo. Most of the bunnies are obvious and easy to see, but some are very cleverly hidden.
Left: June 1989. The bunny logo is disguised as ribbon decoration on the model’s left thigh.
Right: August 1989. The bunny logo is on the newspaper as part of a graph near the business woman’s left thumb.
Demo program. Very long and very complicated. Replace “lt” (less than), “gt”, “lte”, “gte” with Boolean operator symbols (my blog editor chokes on symbols).
using System.IO;
namespace MatrixMoorePenrosePseudoInverseThreeDifferentQR
{
internal class Program
{
static void Main(string[] args)
{
Console.WriteLine("\nBegin relaxed Moore-Penrose" +
" pseudo-inverse using QR decomp Householder ");
//Console.WriteLine("\nBegin relaxed Moore-Penrose" +
// " pseudo-inverse using QR decomp Gram-Schmidt ");
//Console.WriteLine("\nBegin relaxed Moore-Penrose" +
// " pseudo-inverse using QR decomp Givens ");
int minRows = 100; int maxRows = 1000;
int minCols = 2; int maxCols = 20;
Random rnd = new Random(0);
int nTrials = 5000;
int nPass = 4990; int nFail = 0;
for (int i = 4990; i "lt" nTrials; ++i)
{
Console.WriteLine("\n==========");
Console.WriteLine("trial " + i);
double[][] A =
MakeRandomMatrix(minRows, maxRows,
minCols, maxCols, rnd);
double[][] Apinv = QRHouseholder.MatPseudoInv(A);
// double[][] Apinv = QRGramSchmidt.MatPseudoInv(A);
// double[][] Apinv = QRGivens.MatPseudoInv(A);
// check (A * Apinv) * A = A
double[][] AApinv = MatProduct(A, Apinv);
double[][] C = MatProduct(AApinv, A);
//Console.WriteLine("\nA = ");
//MatShow(A, 4, 9);
//Console.WriteLine("\nApinv = ");
//MatShow(Apinv, 4, 9);
//Console.WriteLine("\n(A * Apinv) * A = ");
//MatShow(C, 4, 9);
if (MatAreEqual(C, A, 1.0e-8) == true)
{
Console.WriteLine("pass");
++nPass;
}
else
{
Console.WriteLine("FAIL");
Console.WriteLine("nRows = " + A.Length +
" nCols = " + A[0].Length);
//Console.ReadLine();
++nFail;
}
Console.WriteLine("==========");
// Console.ReadLine();
} // nTrials
Console.WriteLine("\nNumber pass = " + nPass);
Console.WriteLine("Number fail = " + nFail);
Console.WriteLine("\nEnd testing ");
Console.ReadLine();
} // Main
// ------------------------------------------------------
// helpers:
// ------------------------------------------------------
public static double[][] MakeRandomMatrix(int minRows,
int maxRows, int minCols, int maxCols, Random rnd)
{
int nRows = rnd.Next(minRows, maxRows); // [2,maxN)
int nCols = rnd.Next(minCols, maxCols);
double[][] result = new double[nRows][];
for (int i = 0; i "lt" nRows; ++i)
result[i] = new double[nCols];
double lo = -10.0; double hi = 10.0;
for (int r = 0; r "lt" nRows; ++r)
for (int c = 0; c "lt" nCols; ++c)
result[r][c] = (hi - lo) * rnd.NextDouble() + lo;
return result;
}
// ------------------------------------------------------
public static double[][] MatProduct(double[][] A,
double[][] B)
{
int aRows = A.Length; int aCols = A[0].Length;
int bRows = B.Length; int bCols = B[0].Length;
if (aCols != bRows)
throw new Exception("Non-conformable matrices");
double[][] result = new double[aRows][];
for (int i = 0; i "lt" aRows; ++i)
result[i] = new double[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] += A[i][k] * B[k][j];
return result;
}
// ------------------------------------------------------
public static void MatShow(double[][] m, int dec, int wid)
{
int nRows = m.Length; int nCols = m[0].Length;
double small = 1.0 / Math.Pow(10, dec);
for (int i = 0; i "lt" nRows; ++i)
{
for (int j = 0; j "lt" nCols; ++j)
{
double v = m[i][j];
if (Math.Abs(v) "lt" small) v = 0.0;
Console.Write(v.ToString("F" + dec).
PadLeft(wid));
}
Console.WriteLine("");
}
}
// ------------------------------------------------------
public static bool MatAreEqual(double[][] A, double[][] B,
double epsilon)
{
int nRows = A.Length; int nCols = A[0].Length;
for (int i = 0; i "lt" nRows; ++i)
for (int j = 0; j "lt" nCols; ++j)
if (Math.Abs( A[i][j] - B[i][j] ) "gt" epsilon)
return false;
return true;
}
// ------------------------------------------------------
} // Program
// ========================================================
public class QRHouseholder
{
// container class for relaxed Moore-Penrose
// pseudo-inverse using QR decomposition with QR
// Householder algorithm
public static double[][] MatPseudoInv(double[][] M)
{
// relaxed Moore-Penrose pseudo-inverse
// using QR Householder decomposition
// A = Q*R, pinv(A) = inv(R) * trans(Q)
int nr = M.Length; int nc = M[0].Length; // aka m, n
if (nr "lt" nc)
Console.WriteLine("ERROR: Works only m "gte" n");
double[][] Q; double[][] R;
MatDecompQR(M, out Q, out R); // Householder
double[][] Rinv = MatInvUpperTri(R); // std algo
double[][] Qinv = MatTranspose(Q); // is inv(Q)
double[][] result = MatProduct(Rinv, Qinv);
return result;
}
// ------------------------------------------------------
private static void MatDecompQR(double[][] A,
out double[][] Q, out double[][] R)
{
// reduced QR decomp using Householder reflections
int m = A.Length; int n = A[0].Length;
//if (m "lt" n)
// throw new Exception("No rows less than cols");
double[][] QQ = MatIdentity(m); // working Q
double[][] RR = MatCopy(A); // working R
for (int k = 0; k "lt" n; ++k)
{
double[] x = MatColToEnd(RR, k); // RR[k:, k]
double normx = VecNorm(x);
if (Math.Abs(normx) "lt" 1.0e-12) continue;
double[] v = VecCopy(x);
double sign;
if (x[0] "gte" 0.0) sign = 1.0; else sign = -1.0;
v[0] += sign * normx;
double normv = VecNorm(v);
for (int i = 0; i "lt" v.Length; ++i)
v[i] /= normv;
// apply reflection to R
double[][] tmp1 =
MatRowColToEnd(RR, k); // RR[k:, k:]
double[] tmp2 = VecMatProduct(v, tmp1);
double[][] tmp3 = VecOuter(v, tmp2);
double[][] B = MatScalarMult(2.0, tmp3);
MatSubtractCorner(RR, k, B); // R[k:, k:] -= B
// accumulate Q
double[][] tmp4 =
MatExtractAllRowsColToEnd(QQ, k); // QQ[:, k:]
double[] tmp5 = MatVecProduct(tmp4, v);
double[][] tmp6 = VecOuter(tmp5, v);
double[][] C = MatScalarMult(2.0, tmp6);
MatSubtractRows(QQ, k, C); // QQ[:, k:] -= C
} // for k
// return parts of QQ and RR
Q = MatExtractFirstCols(QQ, n); // Q[:, :n]
R = MatExtractFirstRows(RR, n); // R[:n, :]
} // MatDecompQR()
// ------------------------------------------------------
private static double[][] MatInvUpperTri(double[][] A)
{
// used to invert R from QR
int n = A.Length; // must be square matrix
double[][] result = MatIdentity(n);
for (int k = 0; k "lt" n; ++k)
{
for (int j = 0; j "lt" n; ++j)
{
for (int i = 0; i "lt" k; ++i)
{
result[j][k] -= result[j][i] * A[i][k];
}
result[j][k] /= (A[k][k] + 1.0e-8); // avoid 0
}
}
return result;
}
// ------------------------------------------------------
private static double[][] MatTranspose(double[][] M)
{
int nRows = M.Length; int nCols = M[0].Length;
double[][] result = MatMake(nCols, nRows); // note
for (int i = 0; i "lt" nRows; ++i)
for (int j = 0; j "lt" nCols; ++j)
result[j][i] = M[i][j];
return result;
}
// ------------------------------------------------------
private static double[][] MatProduct(double[][] A,
double[][] B)
{
int aRows = A.Length; int aCols = A[0].Length;
int bRows = B.Length; int bCols = B[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] += A[i][k] * B[k][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[][] 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[] MatColToEnd(double[][] A,
int k)
{
// last part of col k, from [k,k] to end
int m = A.Length; int n = A[0].Length;
double[] result = new double[m - k];
for (int i = 0; i "lt" m - k; ++i)
result[i] = A[i + k][k];
return result;
}
// ------------------------------------------------------
private static double[][] MatRowColToEnd(double[][] A,
int k)
{
// block of A, row k to end, col k to end
// A[k:, k:]
int m = A.Length; int n = A[0].Length;
double[][] result = MatMake(m - k, n - k);
for (int i = 0; i "lt" m - k; ++i)
for (int j = 0; j "lt" n - k; ++j)
result[i][j] = A[i + k][j + k];
return result;
}
// ------------------------------------------------------
private static double[] VecCopy(double[] v)
{
double[] result = new double[v.Length];
for (int i = 0; i "lt" v.Length; ++i)
result[i] = v[i];
return result;
}
// ------------------------------------------------------
private static double[] VecMatProduct(double[] v,
double[][] A)
{
// v * A
int m = A.Length; int n = A[0].Length;
double[] result = new double[n];
for (int j = 0; j "lt" n; ++j)
for (int i = 0; i "lt" m; ++i)
result[j] += v[i] * A[i][j];
return result;
}
// ------------------------------------------------------
private static double[][] VecOuter(double[] v1,
double[] v2)
{
// vector outer product
double[][] result = MatMake(v1.Length, v2.Length);
for (int i = 0; i "lt" v1.Length; ++i)
for (int j = 0; j "lt" v2.Length; ++j)
result[i][j] = v1[i] * v2[j];
return result;
}
// ------------------------------------------------------
private static double[][] MatScalarMult(double x,
double[][] A)
{
// x * A element-wise
int m = A.Length; int n = A[0].Length;
double[][] result = MatMake(m, n);
for (int i = 0; i "lt" m; ++i)
for (int j = 0; j "lt" n; ++j)
result[i][j] = x * A[i][j];
return result;
}
// ------------------------------------------------------
private static void MatSubtractCorner(double[][] A,
int k, double[][] C)
{
// subtract elements in C from lower right A
// A[k:, k:] -= C
int m = A.Length; int n = A[0].Length;
for (int i = 0; i "lt" m - k; ++i)
for (int j = 0; j "lt" n - k; ++j)
A[i + k][j + k] -= C[i][j];
return; // modify A in-place
}
// ------------------------------------------------------
private static double[][]
MatExtractAllRowsColToEnd(double[][] A, int k)
{
// A[:, k:]
int m = A.Length; int n = A[0].Length;
double[][] result = MatMake(m, n - k);
for (int i = 0; i "lt" m; ++i)
for (int j = 0; j "lt" n - k; ++j)
result[i][j] = A[i][j + k];
return result;
}
// -----------------------------------------------------
private static double[] MatVecProduct(double[][] A,
double[] v)
{
// A * v
int m = A.Length; int n = A[0].Length;
double[] result = new double[m];
for (int i = 0; i "lt" m; ++i)
for (int j = 0; j "lt" n; ++j)
result[i] += A[i][j] * v[j];
return result;
}
// ------------------------------------------------------
private static void MatSubtractRows(double[][] A,
int k, double[][] C)
{
// A[:, k:] -= C
int m = A.Length; int n = A[0].Length;
for (int i = 0; i "lt" m; ++i)
for (int j = 0; j "lt" n - k; ++j)
A[i][j + k] -= C[i][j];
return; // modify A in-place
}
// ------------------------------------------------------
private static double[][]
MatExtractFirstCols(double[][] A, int k)
{
// A[:, :n]
int m = A.Length; int n = A[0].Length;
double[][] result = MatMake(m, k);
for (int i = 0; i "lt" m; ++i)
for (int j = 0; j "lt" k; ++j)
result[i][j] = A[i][j];
return result;
}
// ------------------------------------------------------
private static double[][]
MatExtractFirstRows(double[][] A, int k)
{
// A[:n, :]
int m = A.Length; int n = A[0].Length;
double[][] result = MatMake(k, n);
for (int i = 0; i "lt" k; ++i)
for (int j = 0; j "lt" n; ++j)
result[i][j] = A[i][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);
}
// ------------------------------------------------------
} // MatDecomQR2
// ========================================================
public class QRGramSchmidt
{
// container class for relaxed Moore-Penrose
// pseudo-inverse using QR decomposition with modified
// Gram-Schmidt algorithm
public static double[][] MatPseudoInv(double[][] M)
{
// relaxed Moore-Penrose pseudo-inverse using QR-GS
// A = Q*R, pinv(A) = inv(R) * trans(Q)
int nr = M.Length; int nc = M[0].Length; // aka m, n
if (nr "lt" nc)
Console.WriteLine("ERROR: Works only m "gte" n");
double[][] Q; double[][] R;
MatDecompQR(M, out Q, out R); // Gram-Schmidt
double[][] Rinv = MatInvUpperTri(R); // std algo
double[][] Qinv = MatTranspose(Q); // is inv(Q)
double[][] result = MatProduct(Rinv, Qinv);
return result;
}
// ------------------------------------------------------
private static void MatDecompQR(double[][] A,
out double[][] Q, out double[][] R)
{
// QR decomposition, modified Gram-Schmidt
// 'reduced' mode (all machine learning scenarios)
int m = A.Length; int n = A[0].Length;
if (m "lt" n)
throw new Exception("m must be gte n ");
double[][] QQ = new double[m][]; // working Q mxn
for (int i = 0; i "lt" m; ++i)
QQ[i] = new double[n];
double[][] RR = new double[n][]; // working R nxn
for (int i = 0; i "lt" n; ++i)
RR[i] = new double[n];
for (int k = 0; k "lt" n; ++k) // main loop each col
{
double[] v = new double[m];
for (int i = 0; i "lt" m; ++i) // col k
v[i] = A[i][k];
for (int j = 0; j "lt" k; ++j) // inner loop
{
double[] colj = new double[QQ.Length];
for (int i = 0; i "lt" colj.Length; ++i)
colj[i] = QQ[i][j];
double vecdot = 0.0;
for (int i = 0; i "lt" colj.Length; ++i)
vecdot += colj[i] * v[i];
RR[j][k] = vecdot;
// v = v - (R[j, k] * Q[:, j])
for (int i = 0; i "lt" v.Length; ++i)
v[i] = v[i] - (RR[j][k] * QQ[i][j]);
} // j
double normv = 0.0;
for (int i = 0; i "lt" v.Length; ++i)
normv += v[i] * v[i];
normv = Math.Sqrt(normv);
RR[k][k] = normv;
// Q[:, k] = v / R[k, k]
for (int i = 0; i "lt" QQ.Length; ++i)
QQ[i][k] = v[i] / (RR[k][k] + 1.0e-12);
} // k
Q = QQ;
R = RR;
}
// ------------------------------------------------------
private static double[][] MatInvUpperTri(double[][] A)
{
// used to invert R from QR
int n = A.Length; // must be square matrix
double[][] result = MatIdentity(n);
for (int k = 0; k "lt" n; ++k)
{
for (int j = 0; j "lt" n; ++j)
{
for (int i = 0; i "lt" k; ++i)
{
result[j][k] -= result[j][i] * A[i][k];
}
result[j][k] /= (A[k][k] + 1.0e-8); // avoid 0
}
}
return result;
}
// ------------------------------------------------------
private static double[][] MatTranspose(double[][] M)
{
int nRows = M.Length; int nCols = M[0].Length;
double[][] result = MatMake(nCols, nRows); // note
for (int i = 0; i "lt" nRows; ++i)
for (int j = 0; j "lt" nCols; ++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[][] A,
double[][] B)
{
int aRows = A.Length; int aCols = A[0].Length;
int bRows = B.Length; int bCols = B[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] += A[i][k] * B[k][j];
return result;
}
} // class QRGramSchmidt
// ========================================================
public class QRGivens
{
// container class for relaxed Moore-Penrose
// pseudo-inverse using QR decomposition with Givens
// rotations algorithm
public static double[][] MatPseudoInv(double[][] M)
{
// relaxed Moore-Penrose pseudo-inverse using QR-GS
// A = Q*R, pinv(A) = inv(R) * trans(Q)
int nr = M.Length; int nc = M[0].Length; // aka m, n
if (nr "lt" nc)
Console.WriteLine("ERROR: Works only m "gte" n");
double[][] Q; double[][] R;
MatDecompQR(M, out Q, out R); // Givens
double[][] Rinv = MatInvUpperTri(R); // std algo
double[][] Qinv = MatTranspose(Q); // is inv(Q)
double[][] result = MatProduct(Rinv, Qinv);
return result;
}
// ------------------------------------------------------
private static void MatDecompQR(double[][] A,
out double[][] Q, out double[][] R)
{
// QR decomposition using Givens rotations
// 'reduced' mode (all machine learning scenarios)
int m = A.Length; int n = A[0].Length;
if (m "lt" n)
throw new Exception("m must be gte n ");
int k = Math.Min(m, n);
double[][] QQ = new double[m][]; // working Q mxm
for (int i = 0; i "lt" m; ++i)
QQ[i] = new double[m];
for (int i = 0; i "lt" m; ++i)
QQ[i][i] = 1.0; // Identity
double[][] RR = new double[m][]; // working R mxn
for (int i = 0; i "lt" m; ++i)
RR[i] = new double[n];
for (int i = 0; i "lt" m; ++i)
for (int j = 0; j "lt" n; ++j)
RR[i][j] = A[i][j]; // copy of A
for (int j = 0; j "lt" k; ++j) // outer loop
{
for (int i = m - 1; i "gt" j; --i) // inner loop
{
double a = RR[j][j];
double b = RR[i][j];
//if (b == 0.0) continue; // risky
if (Math.Abs(b) "lt" 1.0e-12) continue;
double r = Math.Sqrt((a * a) + (b * b));
double c = a / r;
double s = -b / r;
for (int col = j; col "lt" n; ++col)
{
double tmp = (c * RR[j][col]) - (s * RR[i][col]);
RR[i][col] = (s * RR[j][col]) + (c * RR[i][col]);
RR[j][col] = tmp;
}
for (int row = 0; row "lt" m; ++row)
{
double tmp = (c * QQ[row][j]) - (s * QQ[row][i]);
QQ[row][i] = (s * QQ[row][j]) + (c * QQ[row][i]);
QQ[row][j] = tmp;
}
} // inner i loop
} // j loop
// Q_reduced = Q[:, :k] // all rows, col 0 to k
// R_reduced = R[:k, :] // rows 0 to k, all cols
double[][] Qred = new double[m][];
for (int i = 0; i "lt" m; ++i)
Qred[i] = new double[k];
for (int i = 0; i "lt" m; ++i)
for (int j = 0; j "lt" k; ++j)
Qred[i][j] = QQ[i][j];
Q = Qred;
double[][] Rred = new double[k][];
for (int i = 0; i "lt" k; ++i)
Rred[i] = new double[n];
for (int i = 0; i "lt" k; ++i)
for (int j = 0; j "lt" n; ++j)
Rred[i][j] = RR[i][j];
R = Rred;
} // MatDecompQR()
// ------------------------------------------------------
private static double[][] MatInvUpperTri(double[][] A)
{
// used to invert R from QR
int n = A.Length; // must be square matrix
double[][] result = MatIdentity(n);
for (int k = 0; k "lt" n; ++k)
{
for (int j = 0; j "lt" n; ++j)
{
for (int i = 0; i "lt" k; ++i)
{
result[j][k] -= result[j][i] * A[i][k];
}
result[j][k] /= (A[k][k] + 1.0e-8); // avoid 0
}
}
return result;
}
// ------------------------------------------------------
private static double[][] MatTranspose(double[][] M)
{
int nRows = M.Length; int nCols = M[0].Length;
double[][] result = MatMake(nCols, nRows); // note
for (int i = 0; i "lt" nRows; ++i)
for (int j = 0; j "lt" nCols; ++j)
result[j][i] = M[i][j];
return result;
}
// ------------------------------------------------------
private static double[][] MatProduct(double[][] A,
double[][] B)
{
int aRows = A.Length; int aCols = A[0].Length;
int bRows = B.Length; int bCols = B[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] += A[i][k] * B[k][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;
}
// ------------------------------------------------------
} // class QRGivens
// ========================================================
} // ns
.NET Test Automation Recipes
Software Testing
SciPy Programming Succinctly
Keras Succinctly
R Programming
2026 Visual Studio Live
2025 Summer MLADS Conference
2026 DevIntersection Conference
2025 Machine Learning Week
2025 Ai4 Conference
2026 G2E Conference
2026 iSC West Conference
You must be logged in to post a comment.