One rainy morning before work, I figured I’d dust off some of my C# matrix inverse code. If you have a square matrix A, the inverse, Ai, is a matrix where A * Ai = I, where I is the identity matrix. (Also, Ai * A = I). For example, if A =
4.00 7.00 2.00 6.00 0.00 3.00 8.00 1.00 9.00
The inverse of A is:
0.014286 0.290476 -0.100000 0.142857 -0.095238 0.000000 -0.028571 -0.247619 0.200000
And A * Ai = I which is:
1.00 0.00 0.00 0.00 1.00 0.00 0.00 0.00 1.00
Computing the inverse of a matrix is extremely difficult. Wikipedia lists over a dozen different algorithms, and each algorithm has several variations. The number of different algorithms is a direct indication of the difficulty of matrix inversion.
I put together a demo of matrix inverse with C# using five different algorithms: 1.) LUP (Croat version), 2. QR (Householder version), 3.) SVD (Jacobi version), 4.) Newton (Pan version), 5.) Cholesky (Banachiewicz version).
Here’s the output of my demo:
Begin C# matrix inverse demo Source matrix A: 4.00 7.00 2.00 6.00 0.00 3.00 8.00 1.00 9.00 ========== Inverse matrix via LUP-Crout: 0.014286 0.290476 -0.100000 0.142857 -0.095238 0.000000 -0.028571 -0.247619 0.200000 Inverse matrix via QR-Householder: 0.014286 0.290476 -0.100000 0.142857 -0.095238 -0.000000 -0.028571 -0.247619 0.200000 Inverse matrix via SVD-Jacobi: 0.014286 0.290476 -0.100000 0.142857 -0.095238 -0.000000 -0.028571 -0.247619 0.200000 Inverse matrix via Newton-Pan: 0.014286 0.290476 -0.100000 0.142857 -0.095238 -0.000000 -0.028571 -0.247619 0.200000 ========== Postive-Definite matrix M: 1.00 0.20 0.30 0.40 0.20 1.00 0.50 0.60 0.30 0.50 1.00 0.70 0.40 0.60 0.70 1.00 Inverse matrix via Cholesky-Banachiewicz: 1.195815 0.082212 -0.059791 -0.485800 0.082212 1.599402 -0.254111 -0.814649 -0.059791 -0.254111 2.002990 -1.225710 -0.485800 -0.814649 -1.225710 2.541106 verified correct End
Not every square matrix has an inverse. There is no one best matrix inverse technique. Each technique has pros and cons. The Cholesky algorithm applies only to symmetric matrices (actually positive-definite) that occur in many algorithms that use kernel matrices, such as kernel ridge regression.
It would literally take an entire book to explain all of these techniques in detail, so I won’t try. The code is here in case someone out there needs it.
Good fun for me on a rainy morning before work.
I like implementing code from scratch instead of using an external library. Using from-scratch code gives me full visibility into the inner workings of the system.
The girl-in-the-glass-tube is a common science fiction theme. Left: “Amazing Stories”, Nov. 1939. Center: “Frozen Limit”, 1950, by John Fearn. Right: “Planet Stories”, Summer 1948.
Demo code. Very long and very complex. I believe it’s all solid. It’s been tested but not exhaustively. 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;
// C# matrix inverse using five different algorithms
// 1. LUP (Crout)
// 2. QR (Householder)
// 3. SVD (Jacobi)
// 4. Newton (Pan)
// 5. Cholesky (Banachiewicz) (positive definite only)
namespace MatrixInverses
{
internal class MatrixInversesProgram
{
static void Main(string[] args)
{
Console.WriteLine("\nBegin C# matrix inverse demo ");
double[][] A = new double[4][];
A[0] = new double[] { 4.0, 7.0, 1.0, 2.0 };
A[1] = new double[] { 6.0, 0.0, 3.0, 5.0 };
A[2] = new double[] { 8.0, 1.0, 9.0, 2.0 };
A[3] = new double[] { 2.0, 5.0, 6.0, -3.0 };
// to load matrix from file:
// double[][] A = MatLoad("my_data.txt",
// new int[] {0, 1, 2, 3}, ',', "#");
Console.WriteLine("\nSource matrix A: ");
MatShow(A, 2, 6);
Console.WriteLine("\n==========");
double[][] AiLUP = MatInverseLUP(A);
Console.WriteLine("\nInverse matrix via" +
" LUP-Crout: ");
MatShow(AiLUP, 6, 11);
double[][] AiQR = MatInverseQR(A);
Console.WriteLine("\nInverse matrix via" +
" QR-Householder: ");
MatShow(AiQR, 6, 11);
double[][] AiSVD = MatInverseSVD(A);
Console.WriteLine("\nInverse matrix via" +
" SVD-Jacobi: ");
MatShow(AiSVD, 6, 11);
double[][] AiNewton = MatInverseNewton(A);
Console.WriteLine("\nInverse matrix via" +
" Newton-Pan: ");
MatShow(AiNewton, 6, 11);
double[][] M = new double[4][]; // positive definite
M[0] = new double[] { 1.0, 0.2, 0.3, 0.4 };
M[1] = new double[] { 0.2, 1.0, 0.5, 0.6 };
M[2] = new double[] { 0.3, 0.5, 1.0, 0.7 };
M[3] = new double[] { 0.4, 0.6, 0.7, 1.0 };
Console.WriteLine("\n==========");
Console.WriteLine("\nPostive-Definite matrix M: ");
MatShow(M, 2, 6);
double[][] MiC = MatInverseCholesky(M);
Console.WriteLine("\nInverse matrix via" +
" Cholesky-Banachiewicz: ");
MatShow(MiC, 6, 11);
double[][] check = MatProduct(M, MiC);
double[][] I = MatIdentity(M.Length);
if (MatAreEqual(check, I, 1.0e-12) == true)
Console.WriteLine("verified correct ");
Console.WriteLine("\nEnd ");
Console.ReadLine();
} // Main
// ------------------------------------------------------
// helpers for Main():
static double[][] MatIdentity(int n)
{
double[][] result = new double[n][];
for (int i = 0; i "lt" n; ++i)
result[i] = new double[n];
for (int i = 0; i "lt" n; ++i)
result[i][i] = 1.0;
return result;
}
static bool MatAreEqual(double[][] A,
double[][] B, double eps)
{
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 (Math.Abs(A[i][j] - B[i][j]) "gt" eps)
return false;
return true;
}
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("");
}
}
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;
}
static double[][] MatLoad(string fn, int[] usecols,
char sep, string comment)
{
List"lt"double[]"gt" result =
new List"lt"double[]"gt"();
string line = "";
FileStream ifs = new FileStream(fn, FileMode.Open);
StreamReader sr = new StreamReader(ifs);
while ((line = sr.ReadLine()) != null)
{
if (line.StartsWith(comment) == true)
continue;
string[] tokens = line.Split(sep);
List"lt"double"gt" lst = new List"lt"double"gt"();
for (int j = 0; j "lt" usecols.Length; ++j)
lst.Add(double.Parse(tokens[usecols[j]]));
double[] row = lst.ToArray();
result.Add(row);
}
sr.Close(); ifs.Close();
return result.ToArray();
}
// ======================================================
// 1. LUP - Crout
// ======================================================
static double[][] MatInverseLUP(double[][] M)
{
// LUP using Crout's algorithm
// A = P * L * U
// Ainv = inv(U) * inv(L) * inv(P)
// note: inv(P) = P
int n = M.Length;
double[][] lum;
int[] perm;
int sgn = MatDecomposeLUP(M, out lum, out perm);
// check if determinant is 0
double det = (double)sgn;
for (int i = 0; i "lt" n; ++i)
det *= lum[i][i];
if (det == 0.0)
{
Console.WriteLine("Inverse doesn't exist ");
Console.ReadLine();
}
// call two helper functions . .
double[][] L = ExtractLower(lum); // 1s on diag
double[][] U = ExtractUpper(lum); // upper tri
// convert perm vector to P matrix
double[][] P = MatMake(n, n);
for (int i = 0; i "lt" n; ++i)
P[i][perm[i]] = 1.0;
// compute Li and Ui via helpers
double[][] Li = MatInverseLowerTri(L);
double[][] Ui = MatInverseUpperTri(U);
// compute inverse Ai = Ui * Li * Pi
double[][] result = MatProduct(Ui, MatProduct(Li, P));
return result;
// ----------------------------------------------------
// helpers: MatDecomposeLUP, MatMake, ExtractLower,
// ExtractUpper, MatInverseLowerTri,
// MatInverseUpperTri, MatIdentity, MatProduct
// ----------------------------------------------------
static int MatDecomposeLUP(double[][] m,
out double[][] lum, out int[] perm)
{
// Crout's LU decomposition for matrix determinant
// and inverse
// stores combined lower and upper in lum[][]
// stores row permuations into perm[]
// returns +1 or -1 according to even or odd
// number of row permutations
// lower gets dummy 1.0s on diagonal (0.0s above)
// upper gets lum values on diagonal (0.0s below)
int toggle = +1; // even (+1) or odd (-1) perms
int n = m.Length;
// make a copy of m[][] into result lu[][]
lum = MatMake(n, n);
for (int i = 0; i "lt" n; ++i)
for (int j = 0; j "lt" n; ++j)
lum[i][j] = m[i][j];
// make perm[]
perm = new int[n];
for (int i = 0; i "lt" n; ++i)
perm[i] = i;
for (int j = 0; j "lt" n - 1; ++j)
{
double max = Math.Abs(lum[j][j]);
int piv = j;
for (int i = j + 1; i "lt" n; ++i) // pivot index
{
double xij = Math.Abs(lum[i][j]);
if (xij "gt" max)
{
max = xij;
piv = i;
}
} // i
if (piv != j)
{
double[] tmp = lum[piv]; // swap rows j, piv
lum[piv] = lum[j];
lum[j] = tmp;
int t = perm[piv]; // swap perm elements
perm[piv] = perm[j];
perm[j] = t;
toggle = -toggle;
}
double xjj = lum[j][j];
if (xjj != 0.0)
{
for (int i = j + 1; i "lt" n; ++i)
{
double xij = lum[i][j] / xjj;
lum[i][j] = xij;
for (int k = j + 1; k "lt" n; ++k)
lum[i][k] -= xij * lum[j][k];
}
}
} // j
return toggle; // for determinant
} // MatDecomposeHelp
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;
}
static double[][] ExtractLower(double[][] lum)
{
// lower part of an LU Crout's decomposition
// (dummy 1.0s on diagonal, 0.0s above)
int n = lum.Length;
double[][] result = MatMake(n, n);
for (int i = 0; i "lt" n; ++i)
{
for (int j = 0; j "lt" n; ++j)
{
if (i == j)
result[i][j] = 1.0;
else if (i "gt" j)
result[i][j] = lum[i][j];
}
}
return result;
}
static double[][] ExtractUpper(double[][] lum)
{
// upper part of an LU (lu values on diagional
// and above, 0.0s below)
int n = lum.Length;
double[][] result = MatMake(n, n);
for (int i = 0; i "lt" n; ++i)
{
for (int j = 0; j "lt" n; ++j)
{
if (i "lte" j)
result[i][j] = lum[i][j];
}
}
return result;
}
static double[][] MatInverseLowerTri(double[][] lower)
{
// inverse of lower triangular non-fancy version
int n = lower.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[k][j] -= result[i][j] * lower[k][i];
}
result[k][j] /= lower[k][k];
}
}
return result;
}
static double[][] MatInverseUpperTri(double[][] U)
{
int n = U.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] * U[i][k];
}
result[j][k] /= (U[k][k] + 1.0e-8);
}
}
return result;
}
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;
}
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;
}
} // MatInverseLUP
// ======================================================
// 2. QR - Householder
// ======================================================
static double[][] MatInverseQR(double[][] M)
{
// inverse using QR decomp (Householder algorithm)
double[][] Q;
double[][] R;
MatDecompQR(M, out Q, out R, true);
// TODO: check if determinant is zero (no inverse)
double[][] Rinv = MatInverseUpperTri(R); // std algo
double[][] Qtrans = MatTranspose(Q); // is inv(Q)
return MatProduct(Rinv, Qtrans);
// ----------------------------------------------------
// helpers: MatDecomposeQR, MatMake, MatTranspose,
// MatInverseUpperTri, MatIdentity, MatCopy, VecNorm,
// VecDot, VecToMat, MatProduct
// ----------------------------------------------------
static void MatDecompQR(double[][] M,
out double[][] Q, out double[][] R, bool reduced)
{
// QR decomposition, Householder algorithm.
// no helper functions except matrix multiply
int m = M.Length;
int n = M[0].Length;
if (m "lt" n)
throw new Exception("Cannot do rows less than cols");
// double[][] QQ = MatIdentity(m); // working Q
double[][] QQ = new double[m][];
for (int ii = 0; ii "lt" m; ++ii)
QQ[ii] = new double[m];
for (int ii = 0; ii "lt" m; ++ii)
QQ[ii][ii] = 1.0;
//double[][] RR = MatCopy(M); // working R
double[][] RR = new double[m][];
for (int ii = 0; ii "lt" m; ++ii)
RR[ii] = new double[n];
for (int ii = 0; ii "lt" m; ++ii)
for (int jj = 0; jj "lt" n; ++jj)
RR[ii][jj] = M[ii][jj];
int end;
if (m == n) end = n - 1;
else end = n;
for (int i = 0; i "lt" end; ++i)
{
// double[][] H = MatIdentity(m);
double[][] H = new double[m][];
for (int ii = 0; ii "lt" m; ++ii)
H[ii] = new double[m];
for (int ii = 0; ii "lt" m; ++ii)
H[ii][ii] = 1.0;
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);
double sum = 0.0;
for (int ii = 0; ii "lt" a.Length; ++ii)
sum += a[ii] * a[ii];
double normA = Math.Sqrt(sum);
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[][] h = new double[a.Length][];
for (int ii = 0; ii "lt" a.Length; ++ii)
h[ii] = new double[a.Length];
for (int ii = 0; ii "lt" a.Length; ++ii)
h[ii][ii] = 1.0;
// double vvDot = VecDot(v, v);
double vvDot = 0.0;
for (int ii = 0; ii "lt" v.Length; ++ii)
vvDot += v[ii] * v[ii];
// double[][] A = VecToMat(v, v.Length, 1
double[][] A = new double[v.Length][];
int ka = 0;
for (int ii = 0; ii "lt" v.Length; ++ii)
A[ii] = new double[1];
for (int ii = 0; ii "lt" A.Length; ++ii)
for (int jj = 0; jj "lt" A[0].Length; ++jj)
A[ii][jj] = v[ka++];
// double[][] B = VecToMat(v, 1, v.Length);
double[][] B = new double[1][];
for (int ii = 0; ii "lt" 1; ++ii)
B[ii] = new double[v.Length];
int kb = 0;
for (int ii = 0; ii "lt" B.Length; ++ii)
for (int jj = 0; jj "lt" B[0].Length; ++jj)
B[ii][jj] = v[kb++];
double[][] AB = MatMultInternal(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 = MatMultInternal(QQ, H);
RR = MatMultInternal(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); // always n
// double[][] Rsquared = MatMake(dim, dim);
double[][] Rsquared = new double[dim][];
for (int ii = 0; ii "lt" dim; ++ii)
Rsquared[ii] = new double[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);
double[][] Qtrimmed = new double[qRows][];
for (int ii = 0; ii "lt" qRows; ++ii)
Qtrimmed[ii] = new double[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;
}
// --------------------------------------------------
// nested helper MatMultInternal
// --------------------------------------------------
static double[][] MatMultInternal(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);
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;
}
} // MatDecompQR()
static double[][] MatInverseUpperTri(double[][] U)
{
int n = U.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] * U[i][k];
}
result[j][k] /= (U[k][k] + 1.0e-8); // no 0
}
}
return result;
}
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;
}
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;
}
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;
}
} // MatInverseQR
// ======================================================
// 3. SVD - Jacobi
// ======================================================
static double[][] MatInverseSVD(double[][] M)
{
// SVD Jacobi algorithm
// A = U * S * Vh
// inv(A) = tr(Vh) * inv(S) * tr(U)
double[][] U; double[][] Vh; double[] s;
MatDecomposeSVD(M, out U, out Vh, out s);
// TODO: check if determinant is zero (no inverse)
// convert s vector to Sinv matrix
int n = M.Length;
double[][] Sinv = MatMake(n, n);
for (int i = 0; i "lt" n; ++i)
Sinv[i][i] = 1.0 / s[i]; // check s[i] != 0
double[][] V = MatTranspose(Vh);
double[][] Uh = MatTranspose(U);
double[][] result = MatProduct(V,
MatProduct(Sinv, Uh));
return result;
// ----------------------------------------------------
static void MatDecomposeSVD(double[][] mat,
out double[][] U, out double[][] Vh,
out double[] s)
{
// assumes source matrix is square
// see github.com/ampl/gsl/blob/master/linalg/svd.c
// "full_matrices=False" version
double EPSILON = 1.0e-15;
double[][] A = MatCopy(mat); // working U
int m = A.Length;
int n = A[0].Length; // check m == n
double[][] Q = MatIdentity(n); // working V
double[] t = new double[n]; // working s
int ct = 1; // rotation counter
int pass = 0;
double tol = 10 * n * EPSILON; // heuristic
int passMax = 5 * n;
if (passMax "lt" 15) passMax = 15; // heuristic
// save the column error estimates
for (int j = 0; j "lt" n; ++j)
{
double[] cj = MatGetColumn(A, j);
double sj = VecNorm(cj);
t[j] = EPSILON * sj;
}
while (ct "gt" 0 && pass "lte" passMax)
{
ct = n * (n - 1) / 2; // rotation counter
for (int j = 0; j "lt" n - 1; ++j)
{
for (int k = j + 1; k "lt" n; ++k)
{
double sin; double cos;
double[] cj = MatGetColumn(A, j);
double[] ck = MatGetColumn(A, k);
double p = 2.0 * VecDot(cj, ck);
double a = VecNorm(cj);
double b = VecNorm(ck);
double q = a * a - b * b;
double v = Hypot(p, q);
double errorA = t[j];
double errorB = t[k];
bool sorted = false;
if (a "gte" b) sorted = true;
bool orthog = false;
if (Math.Abs(p) "lte"
tol * (a * b)) orthog = true;
bool badA = false;
if (a "lt" errorA) badA = true;
bool badB = false;
if (b "lt" errorB) badB = true;
if (sorted == true && (orthog == true ||
badA == true || badB == true))
{
--ct;
continue;
}
// compute rotation angles
if (v == 0 || sorted == false)
{
cos = 0.0;
sin = 1.0;
}
else
{
cos = Math.Sqrt((v + q) / (2.0 * v));
sin = p / (2.0 * v * cos);
}
// apply rotation to A (U)
for (int i = 0; i "lt" m; ++i)
{
double Aik = A[i][k];
double Aij = A[i][j];
A[i][j] = Aij * cos + Aik * sin;
A[i][k] = -Aij * sin + Aik * cos;
}
// update singular values
t[j] = Math.Abs(cos) * errorA +
Math.Abs(sin) * errorB;
t[k] = Math.Abs(sin) * errorA +
Math.Abs(cos) * errorB;
// apply rotation to Q (V)
for (int i = 0; i "lt" n; ++i)
{
double Qij = Q[i][j];
double Qik = Q[i][k];
Q[i][j] = Qij * cos + Qik * sin;
Q[i][k] = -Qij * sin + Qik * cos;
} // i
} // k
} // j
++pass;
} // while
// compute singular values
double prevNorm = -1.0;
for (int j = 0; j "lt" n; ++j)
{
double[] column = MatGetColumn(A, j);
double norm = VecNorm(column);
// check if singular value is zero
if (norm == 0.0 || prevNorm == 0.0
|| (j "gt" 0 && norm "lte" tol * prevNorm))
{
t[j] = 0.0;
for (int i = 0; i "lt" m; ++i)
A[i][j] = 0.0;
prevNorm = 0.0;
}
else
{
t[j] = norm;
for (int i = 0; i "lt" m; ++i)
A[i][j] = A[i][j] * 1.0 / norm;
prevNorm = norm;
}
}
if (ct "gt" 0)
{
Console.WriteLine("Jacobi iterations did not" +
" converge");
}
U = A;
Vh = MatTranspose(Q);
s = t;
// assume M is square so no need
// to sync with default np.linalg.svd() shapes:
// if (m "lt" n)
// {
// U = MatExtractFirstColumns(U, m);
// s = VecExtractFirst(s, m);
// Vh = MatExtractFirstRows(Vh, m);
// }
} // MatDecomposeSVD
// helpers for MatDecomposeSVD():
//
// MatMake(), MatCopy(), MatIdentity(),
// MatGetColumn(), MatExtractFirstColumns(),
// MatExtractFirstRows(), MatTranspose(),
// MatProduct(), VecNorm(), VecDot(),
// Hypot(), VecExtractFirst()
//
static double[][] MatMake(int nr, int nc)
{
double[][] result = new double[nr][];
for (int i = 0; i "lt" nr; ++i)
result[i] = new double[nc];
return result;
}
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;
}
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;
}
static double[] MatGetColumn(double[][] M, int j)
{
int nRows = M.Length;
double[] result = new double[nRows];
for (int i = 0; i "lt" nRows; ++i)
result[i] = M[i][j];
return result;
}
static double[][] MatExtractFirstColumns(double[][] M,
int n)
{
int nRows = M.Length;
// int nCols = src[0].Length;
double[][] result = MatMake(nRows, n);
for (int i = 0; i "lt" nRows; ++i)
for (int j = 0; j "lt" n; ++j)
result[i][j] = M[i][j];
return result;
}
static double[][] MatExtractFirstRows(double[][] M,
int n)
{
// int nRows = src.Length;
int nCols = M[0].Length;
double[][] result = MatMake(n, nCols);
for (int i = 0; i "lt" n; ++i)
for (int j = 0; j "lt" nCols; ++j)
result[i][j] = M[i][j];
return result;
}
static double[][] MatTranspose(double[][] M)
{
int nr = M.Length;
int nc = M[0].Length;
double[][] result = MatMake(nc, nr);
for (int i = 0; i "lt" nr; ++i)
for (int j = 0; j "lt" nc; ++j)
result[j][i] = M[i][j];
return result;
}
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)
for (int j = 0; j "lt" bCols; ++j)
for (int k = 0; k "lt" aCols; ++k)
result[i][j] += A[i][k] * B[k][j];
return result;
}
static double VecNorm(double[] vec)
{
double sum = 0.0;
int n = vec.Length;
for (int i = 0; i "lt" n; ++i)
sum += vec[i] * vec[i];
return Math.Sqrt(sum);
}
static double VecDot(double[] v1, double[] v2)
{
int n = v1.Length;
double sum = 0.0;
for (int i = 0; i "lt" n; ++i)
sum += v1[i] * v2[i];
return sum;
}
static double Hypot(double x, double y)
{
// fancy sqrt(x^2 + y^2)
double xabs = Math.Abs(x);
double yabs = Math.Abs(y);
double min, max;
if (xabs "lt" yabs)
{
min = xabs; max = yabs;
}
else
{
min = yabs; max = xabs;
}
if (min == 0)
return max;
double u = min / max;
return max * Math.Sqrt(1 + u * u);
}
static double[] VecExtractFirst(double[] vec, int n)
{
double[] result = new double[n];
for (int i = 0; i "lt" n; ++i)
result[i] = vec[i];
return result;
}
} // MatInverseSVD
// ======================================================
// 4. Newton - Pan
// ======================================================
static double[][] MatInverseNewton(double[][] M,
int maxIter = 1000, double epsilon = 1.0e-8)
{
// Newton's method with Pan initialization
// X_k+1 = X_k * (2I - A*X_k)
int n = M.Length; // must be square martix
double[][] Xprev = MatStart(M); // Pan algorithm
double[][] Xnew = MatMake(n, n, 0.0);
double[][] I = MatIdentity(n);
double[][] I2 = MatScalarMult(I, 2.0);
int iter = 0;
while (iter "lt" maxIter)
{
Xnew =
MatProduct(Xprev, MatSubtract(I2,
MatProduct(M, Xprev)));
Xprev = MatCopyOf(Xnew);
if (iter % 10 == 0)
{
double[][] check = MatProduct(M, Xnew); //
if (MatAreEqual(check, I, epsilon) == true)
return Xnew;
}
++iter;
} // while
Console.WriteLine("WARN: no converge ");
return Xnew;
// ----------------------------------------------------
// nested helpers: MatMake(), MatStart(), MatCopyOf(),
// MatTranspose(), MatAreEqual(), MatProduct(),
// MatSubtract(), MatIdentity(), MatScalarMult()
// ----------------------------------------------------
static double[][] MatMake(int nRows, int nCols,
double v)
{
double[][] result = new double[nRows][];
for (int i = 0; i "lt" nRows; ++i)
result[i] = new double[nCols];
for (int i = 0; i "lt" nRows; ++i)
for (int j = 0; j "lt" nCols; ++j)
result[i][j] = v;
return result;
}
// ----------------------------------------------------
static double[][] MatStart(double[][] m)
{
int n = m.Length;
double maxRowSum = 0.0;
double maxColSum = 0.0;
for (int i = 0; i "lt" n; ++i)
{
double rowSum = 0.0;
for (int j = 0; j "lt" n; ++j)
rowSum += Math.Abs(m[i][j]);
if (rowSum "gt" maxRowSum)
maxRowSum = rowSum;
}
for (int j = 0; j "lt" n; ++j)
{
double colSum = 0.0;
for (int i = 0; i "lt" n; ++i)
colSum += Math.Abs(m[i][j]);
if (colSum "gt" maxColSum)
maxColSum = colSum;
}
double[][] result = MatTranspose(m);
double t = 1.0 / (maxRowSum * maxRowSum + 1.0e-8);
for (int i = 0; i "lt" m.Length; ++i)
for (int j = 0; j "lt" m.Length; ++j)
result[i][j] *= t;
return result;
}
// ----------------------------------------------------
static double[][] MatTranspose(double[][] m)
{
int nr = m.Length; int nc = m[0].Length;
double[][] result = MatMake(nc, nr, 0.0); // 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;
}
// ------------------------------------------------------
static double[][] MatCopyOf(double[][] m)
{
int nRows = m.Length; int nCols = m[0].Length;
double[][] result = MatMake(nRows, nCols, 0.0);
for (int i = 0; i "lt" nRows; ++i)
for (int j = 0; j "lt" nCols; ++j)
result[i][j] = m[i][j];
return result;
}
// ------------------------------------------------------
static bool MatAreEqual(double[][] matA,
double[][] matB, double eps)
{
int nr = matA.Length; int nc = matB[0].Length;
for (int i = 0; i "lt" nr; ++i)
for (int j = 0; j "lt" nc; ++j)
if (Math.Abs(matA[i][j] - matB[i][j]) "gt" eps)
return false;
return true;
}
// ----------------------------------------------------
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, 0.0);
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;
}
// ------------------------------------------------------
static double[][] MatSubtract(double[][] matA,
double[][] matB)
{
// matA - matB
int nRows = matA.Length; int nCols = matA[0].Length;
double[][] result = MatMake(nRows, nCols, 0.0);
for (int i = 0; i "lt" nRows; ++i)
for (int j = 0; j "lt" nCols; ++j)
result[i][j] = matA[i][j] - matB[i][j];
return result;
}
// ----------------------------------------------------
static double[][] MatIdentity(int n)
{
double[][] result = MatMake(n, n, 0.0);
for (int i = 0; i "lt" n; ++i)
result[i][i] = 1.0;
return result;
}
// ----------------------------------------------------
static double[][] MatScalarMult(double[][] m, double u)
{
int nRows = m.Length; int nCols = m[0].Length;
double[][] result = MatMake(nRows, nCols, 0.0);
for (int i = 0; i "lt" nRows; ++i)
for (int j = 0; j "lt" nCols; ++j)
result[i][j] = u * m[i][j];
return result;
}
} // MatInverseNewton()
// ======================================================
// 5. Cholesky - Banachiewicz
// ======================================================
static double[][] MatInverseCholesky(double[][] M)
{
double[][] L = MatDecompCholesky(M);
double[][] result = MatInverseFromCholesky(L);
return result;
// ----------------------------------------------------
// helpers: MatDecompCholesky,
// MatInverseFromCholesky, MatIdentity, MatMake
// ----------------------------------------------------
static double[][] MatDecompCholesky(double[][] m)
{
// Cholesky decomposition (Banachiewicz algorithm)
// m is square, symmetric, positive definite
int n = m.Length;
double[][] result = MatMake(n, n); // all 0.0
for (int i = 0; i "lt" n; ++i)
{
for (int j = 0; j "lte" i; ++j)
{
double sum = 0.0;
for (int k = 0; k "lt" j; ++k)
sum += result[i][k] * result[j][k];
if (i == j)
{
double tmp = m[i][i] - sum;
if (tmp "lt" 0.0)
throw new
Exception("MatDecompCholesky fatal");
result[i][j] = Math.Sqrt(tmp);
}
else
{
if (result[j][j] == 0.0)
throw new
Exception("MatDecompCholesky fatal ");
result[i][j] =
(1.0 / result[j][j] * (m[i][j] - sum));
}
} // j
} // i
return result;
} // MatDecompCholesky
// ----------------------------------------------------
static double[][] MatInverseFromCholesky(double[][] L)
{
// L is a lower triangular result of Cholesky decomp
// direct version
int n = L.Length;
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[k][j] -= result[i][j] * L[k][i];
}
result[k][j] /= L[k][k];
}
}
for (int k = n - 1; k "gte" 0; --k)
{
for (int j = 0; j "lt" n; j++)
{
for (int i = k + 1; i "lt" n; i++)
{
result[k][j] -= result[i][j] * L[i][k];
}
result[k][j] /= L[k][k];
}
}
return result;
}
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;
}
static double[][] MatMake(int rows, int cols)
{
double[][] result = new double[rows][];
for (int i = 0; i "lt" rows; ++i)
result[i] = new double[cols];
return result;
}
} // MatInverseCholesky
// ======================================================
} // Program
} // ns
.NET Test Automation Recipes
Software Testing
SciPy Programming Succinctly
Keras Succinctly
R Programming
2026 Visual Studio Live
2025 Summer MLADS Conference
2025 DevIntersection Conference
2025 Machine Learning Week
2025 Ai4 Conference
2025 G2E Conference
2025 iSC West Conference
You must be logged in to post a comment.