feat: add column pivoting to the QR decomposition

README.md

@@ -247,6 +247,21 @@ var Q = QR.orthogonalMatrix;
 var R = QR.upperTriangularMatrix;
 // So you have the QR decomposition. If you multiply Q by R, you'll see that A = Q.R, with Q orthogonal and R upper triangular
 ```
+Pass `pivoting` to move the column carrying the largest remaining norm into place at each step. The diagonal of `R` then comes out non increasing, `rank` becomes reliable, and a least square problem whose matrix is not of full rank can be solved.
+```js
+var A = new Matrix([
+  [1, 2, 3],
+  [2, 4, 6],
+  [1, 1, 1],
+  [3, 5, 7],
+]); // the third column is 2 * second - first, so the rank is 2
+
+var QR = new QrDecomposition(A, { pivoting: true });
+var rank = QR.rank;                     // rank = 2
+var permutation = QR.columnPermutationVector; // the column of A sitting at each position, so A[:, permutation] = Q.R
+var x = QR.solve(Matrix.columnVector([1, 2, 3, 4]));
+// x holds at most `rank` non zero components, the ones the pivoting dropped are pinned to zero
+```
 ##### LU Decomposition
 ```js
 var A = new Matrix([

matrix.d.ts

@@ -1584,18 +1584,45 @@ export class LuDecomposition {
 
 export { LuDecomposition as LU };
 
+export interface IQROptions {
+  /**
+   * Move the column carrying the largest remaining norm into place at each step, which
+   * makes the diagonal of R non increasing and reveals the rank of the input.
+   * Needed to solve a least squares problem whose matrix is not of full rank.
+   * Requires a matrix with at least as many rows as columns.
+   * @default `false`
+   */
+  pivoting?: boolean;
+}
+
 /**
  * @link https://github.com/lutzroeder/Mapack/blob/master/Source/QrDecomposition.cs
  */
 export class QrDecomposition {
-  constructor(value: MaybeMatrix);
+  constructor(value: MaybeMatrix, options?: IQROptions);
   isFullRank(): boolean;
 
   /**
-   * Solve a problem of least square (Ax=b) by using the QR decomposition. Useful when A is rectangular, but not working when A is singular.
+   * The number of diagonal entries of R that stay above a relative tolerance.
+   * Reliable when the decomposition was built with `pivoting`, since the pivoting is
+   * what pushes the negligible entries to the end of the diagonal.
+   */
+  readonly rank: number;
+
+  /**
+   * The column of the input sitting at each position of the decomposition, so that
+   * `A[:, columnPermutationVector] = Q R`. The identity without `pivoting`.
+   */
+  readonly columnPermutationVector: number[];
+
+  /**
+   * Solve a problem of least square (Ax=b) by using the QR decomposition. Useful when A is rectangular.
    * Example : We search to approximate x, with A matrix shape m*n, x vector size n, b vector size m (m > n). We will use :
    * var qr = QrDecomposition(A);
    * var x = qr.solve(b);
+   * Without `pivoting` a matrix that is not of full rank is refused.
+   * With `pivoting` such a matrix gives the basic solution, the one holding at most `rank`
+   * non zero components, the rest pinned to zero.
    * @param value - Matrix 1D which is the vector b (in the equation Ax = b).
    * @returns - The vector x.
    */

src/__tests__/decompositions/qrPivoting.test.js

@@ -0,0 +1,169 @@
+import { describe, it, expect } from 'vitest';
+
+import { Matrix, QrDecomposition } from '../..';
+
+// the third column is 2 * second - first, so the rank is 2
+const rankDeficient = new Matrix([
+  [1, 2, 3],
+  [2, 4, 6],
+  [1, 1, 1],
+  [3, 5, 7],
+]);
+const fullRank = new Matrix([
+  [1, 2],
+  [3, 4],
+  [5, 7],
+]);
+
+function maxAbsDifference(a, b) {
+  let largest = 0;
+  for (let i = 0; i < a.rows; i++) {
+    for (let j = 0; j < a.columns; j++) {
+      largest = Math.max(largest, Math.abs(a.get(i, j) - b.get(i, j)));
+    }
+  }
+  return largest;
+}
+
+function permuteColumns(matrix, permutation) {
+  const result = new Matrix(matrix.rows, matrix.columns);
+  for (let j = 0; j < permutation.length; j++) {
+    result.setColumn(j, matrix.getColumn(permutation[j]));
+  }
+  return result;
+}
+
+describe('QR with column pivoting', () => {
+  it('factors the permuted input', () => {
+    const qr = new QrDecomposition(rankDeficient, { pivoting: true });
+    const permuted = permuteColumns(rankDeficient, qr.columnPermutationVector);
+    const product = qr.orthogonalMatrix.mmul(qr.upperTriangularMatrix);
+    expect(maxAbsDifference(permuted, product)).toBeLessThan(1e-12);
+  });
+
+  it('orders the diagonal of R by decreasing magnitude', () => {
+    const qr = new QrDecomposition(rankDeficient, { pivoting: true });
+    const diagonal = qr.upperTriangularMatrix.diag().map(Math.abs);
+    for (let i = 1; i < diagonal.length; i++) {
+      expect(diagonal[i]).toBeLessThanOrEqual(diagonal[i - 1] + 1e-12);
+    }
+  });
+
+  it('reports the rank of a deficient matrix', () => {
+    const qr = new QrDecomposition(rankDeficient, { pivoting: true });
+    expect(qr.rank).toBe(2);
+    expect(qr.isFullRank()).toBe(false);
+  });
+
+  it('reports the rank of a full rank matrix', () => {
+    const qr = new QrDecomposition(fullRank, { pivoting: true });
+    expect(qr.rank).toBe(2);
+    expect(qr.isFullRank()).toBe(true);
+  });
+
+  it('reports a rank of zero for a matrix of zeros', () => {
+    const qr = new QrDecomposition(new Matrix(3, 2), { pivoting: true });
+    expect(qr.rank).toBe(0);
+  });
+
+  it('returns a permutation of the column indices', () => {
+    const permutation = new QrDecomposition(rankDeficient, {
+      pivoting: true,
+    }).columnPermutationVector;
+    expect([...permutation].sort()).toStrictEqual([0, 1, 2]);
+  });
+
+  it('leaves the permutation as the identity without pivoting', () => {
+    expect(
+      new QrDecomposition(rankDeficient).columnPermutationVector,
+    ).toStrictEqual([0, 1, 2]);
+  });
+
+  it('hands back a copy of the permutation', () => {
+    const qr = new QrDecomposition(rankDeficient, { pivoting: true });
+    const first = qr.columnPermutationVector;
+    first[0] = 99;
+    expect(qr.columnPermutationVector[0]).not.toBe(99);
+  });
+
+  it('rejects a wide matrix instead of indexing past its rows', () => {
+    expect(
+      () =>
+        new QrDecomposition(
+          [
+            [1, 2, 3],
+            [4, 5, 6],
+          ],
+          { pivoting: true },
+        ),
+    ).toThrow(/^Matrix must have at least as many rows as columns$/);
+  });
+});
+
+describe('QR with column pivoting solves a rank deficient least squares', () => {
+  const b = Matrix.columnVector([1, 2, 3, 4]);
+
+  it('satisfies the normal equations', () => {
+    const x = new QrDecomposition(rankDeficient, { pivoting: true }).solve(b);
+    // a least squares solution leaves a residual orthogonal to every column
+    const residual = rankDeficient.mmul(x).sub(b);
+    const projected = rankDeficient.transpose().mmul(residual);
+    expect(maxAbsDifference(projected, new Matrix(3, 1))).toBeLessThan(1e-10);
+  });
+
+  it('pins the dropped components to zero', () => {
+    const qr = new QrDecomposition(rankDeficient, { pivoting: true });
+    const x = qr.solve(b);
+    const permutation = qr.columnPermutationVector;
+    for (let k = qr.rank; k < permutation.length; k++) {
+      expect(x.get(permutation[k], 0)).toBe(0);
+    }
+  });
+
+  it('returns one row per column of the input', () => {
+    const x = new QrDecomposition(rankDeficient, { pivoting: true }).solve(b);
+    expect(x.rows).toBe(3);
+    expect(x.columns).toBe(1);
+  });
+
+  it('solves several right hand sides at once', () => {
+    const rhs = new Matrix([
+      [1, 4],
+      [2, 3],
+      [3, 2],
+      [4, 1],
+    ]);
+    const x = new QrDecomposition(rankDeficient, { pivoting: true }).solve(rhs);
+    expect(x.rows).toBe(3);
+    expect(x.columns).toBe(2);
+    const residual = rankDeficient.mmul(x).sub(rhs);
+    const projected = rankDeficient.transpose().mmul(residual);
+    expect(maxAbsDifference(projected, new Matrix(3, 2))).toBeLessThan(1e-10);
+  });
+
+  it('still refuses a deficient matrix without pivoting', () => {
+    const singular = new Matrix([
+      [1, 1],
+      [0, 0],
+    ]);
+    expect(() =>
+      new QrDecomposition(singular).solve(Matrix.columnVector([1, 0])),
+    ).toThrow(/^Matrix is rank deficient$/);
+  });
+});
+
+describe('QR with column pivoting on a full rank input', () => {
+  const b = Matrix.columnVector([1, 2, 3]);
+
+  it('agrees with the unpivoted solution', () => {
+    const plain = new QrDecomposition(fullRank).solve(b);
+    const pivoted = new QrDecomposition(fullRank, { pivoting: true }).solve(b);
+    expect(maxAbsDifference(plain, pivoted)).toBeLessThan(1e-10);
+  });
+
+  it('rejects a pivoting option that is not a boolean', () => {
+    expect(() => new QrDecomposition(fullRank, { pivoting: 1 })).toThrow(
+      /^pivoting must be a boolean$/,
+    );
+  });
+});

src/dc/qr.js

@@ -4,16 +4,53 @@ import WrapperMatrix2D from '../wrap/WrapperMatrix2D';
 import { hypotenuse } from './util';
 
 export default class QrDecomposition {
-  constructor(value) {
+  constructor(value, options = {}) {
+    const { pivoting = false } = options;
+    if (typeof pivoting !== 'boolean') {
+      throw new TypeError('pivoting must be a boolean');
+    }
     value = WrapperMatrix2D.checkMatrix(value);
+    if (pivoting && value.rows < value.columns) {
+      throw new RangeError('Matrix must have at least as many rows as columns');
+    }
 
     let qr = value.clone();
     let m = value.rows;
     let n = value.columns;
     let rdiag = new Float64Array(n);
     let i, j, k, s;
 
+    // columnPermutation[j] holds the column of the input sitting at position j
+    let columnPermutation = new Array(n);
+    for (k = 0; k < n; k++) {
+      columnPermutation[k] = k;
+    }
+
     for (k = 0; k < n; k++) {
+      if (pivoting) {
+        // move the column with the largest remaining norm into position k, so
+        // that the diagonal of R comes out non increasing and a rank deficient
+        // input pushes its negligible values to the end
+        let best = k;
+        let bestNorm = -1;
+        for (j = k; j < n; j++) {
+          let candidate = 0;
+          for (i = k; i < m; i++) {
+            candidate = hypotenuse(candidate, qr.get(i, j));
+          }
+          if (candidate > bestNorm) {
+            bestNorm = candidate;
+            best = j;
+          }
+        }
+        if (best !== k) {
+          qr.swapColumns(k, best);
+          const swap = columnPermutation[k];
+          columnPermutation[k] = columnPermutation[best];
+          columnPermutation[best] = swap;
+        }
+      }
+
       let nrm = 0;
       for (i = k; i < m; i++) {
         nrm = hypotenuse(nrm, qr.get(i, k));
@@ -42,6 +79,27 @@ export default class QrDecomposition {
 
     this.QR = qr;
     this.Rdiag = rdiag;
+    this.pivoting = pivoting;
+    this.columnPermutation = columnPermutation;
+  }
+
+  get columnPermutationVector() {
+    return this.columnPermutation.slice();
+  }
+
+  get rank() {
+    const n = this.QR.columns;
+    let largest = 0;
+    for (let i = 0; i < n; i++) {
+      largest = Math.max(largest, Math.abs(this.Rdiag[i]));
+    }
+    if (largest === 0) return 0;
+    const tolerance = Math.max(this.QR.rows, n) * largest * Number.EPSILON;
+    let rank = 0;
+    for (let i = 0; i < n; i++) {
+      if (Math.abs(this.Rdiag[i]) > tolerance) rank++;
+    }
+    return rank;
   }
 
   solve(value) {
@@ -53,15 +111,19 @@ export default class QrDecomposition {
     if (value.rows !== m) {
       throw new Error('Matrix row dimensions must agree');
     }
-    if (!this.isFullRank()) {
-      throw new Error('Matrix is rank deficient');
-    }
 
     let count = value.columns;
     let X = value.clone();
     let n = qr.columns;
     let i, j, k, s;
 
+    // without pivoting a rank deficient input has nowhere to put its negligible
+    // columns, so the old refusal stands
+    const rank = this.pivoting ? this.rank : n;
+    if (!this.pivoting && !this.isFullRank()) {
+      throw new Error('Matrix is rank deficient');
+    }
+
     for (k = 0; k < n; k++) {
       for (j = 0; j < count; j++) {
         s = 0;
@@ -74,7 +136,14 @@ export default class QrDecomposition {
         }
       }
     }
-    for (k = n - 1; k >= 0; k--) {
+    // the trailing columns are dropped, which gives the basic solution: the
+    // components they carry are pinned to zero
+    for (k = rank; k < n; k++) {
+      for (j = 0; j < count; j++) {
+        X.set(k, j, 0);
+      }
+    }
+    for (k = rank - 1; k >= 0; k--) {
       for (j = 0; j < count; j++) {
         X.set(k, j, X.get(k, j) / this.Rdiag[k]);
       }
@@ -85,11 +154,25 @@ export default class QrDecomposition {
       }
     }
 
-    return X.subMatrix(0, n - 1, 0, count - 1);
+    const solution = X.subMatrix(0, n - 1, 0, count - 1);
+    if (!this.pivoting) {
+      return solution;
+    }
+    // undo the column swaps so the rows line up with the input columns again
+    const unpermuted = new Matrix(n, count);
+    for (k = 0; k < n; k++) {
+      for (j = 0; j < count; j++) {
+        unpermuted.set(this.columnPermutation[k], j, solution.get(k, j));
+      }
+    }
+    return unpermuted;
   }
 
   isFullRank() {
     let columns = this.QR.columns;
+    if (this.pivoting) {
+      return this.rank === columns;
+    }
     for (let i = 0; i < columns; i++) {
       if (this.Rdiag[i] === 0) {
         return false;