Add Moore-Penrose pseudo inverse and fast non-negative least squares (#455)

* Add Moore-Penrose pseudo inverse and fast non-negative least squares

As we look into to Kubelka-Munk, it became apparent we need some ways to
solve least squares.

Moore-Penrose currently just does basic inverse of tall and wide, but
real tricky matrices will likely not invert and would require a much
more advanced approach. Hopefully we won't need such an approach.

pinv can produce negative values, so if you need a non-negative least
squares, fnnls will do the trick though it is more costly.

Author: facelessuser

coloraide/algebra.py

@@ -159,11 +159,11 @@ def clamp(
         return value
 
 
-def zdiv(a: float, b: float) -> float:
+def zdiv(a: float, b: float, default: float = 0.0) -> float:
     """Protect against zero divide."""
 
     if b == 0:
-        return 0.0
+        return default
     return a / b
 
 
@@ -3592,6 +3592,29 @@ def inv(matrix: MatrixLike | TensorLike) -> Matrix | Tensor:
     return _back_sub_matrix(u, _forward_sub_matrix(l, p, s2), s2)
 
 
+def pinv(a: MatrixLike) -> Matrix:
+    """
+    Compute the (Moore-Penrose) pseudo-inverse of a matrix.
+
+    We currently 'assume' the matrix if full rank. If not, a singular matrix error
+    will be thrown. Such matrices may still be invertible, but they would require
+    a more advanced approach that we do not currently implement.
+
+    Negative results can be returned, use `fnnls` for a non-negative solution (if possible).
+    """
+
+    s = shape(a)
+    if len(s) != 2:
+        raise ValueError('Inputs can only be matrices, vectors or tensors are not allowed')
+
+    t = transpose(a)
+    if s[0] >= s[1]:
+        p = matmul(inv(matmul(t, a, dims=D2)), t, dims=D2)
+    else:
+        p = matmul(t, inv(matmul(a, t, dims=D2)), dims=D2)
+    return p
+
+
 @overload
 def vstack(arrays: Sequence[float | Vector | Matrix]) -> Matrix:
     ...
@@ -3866,3 +3889,88 @@ def inner(a: float | ArrayLike, b: float | ArrayLike) -> float | Array:
 
     # Shape the data.
     return reshape(m, new_shape)  # type: ignore[no-any-return]
+
+
+def fnnls(
+    A: MatrixLike,
+    b: VectorLike,
+    epsilon: float = 1e-12,
+    max_iters: int = 0
+) -> tuple[Vector, float]:
+    """
+    Fast non-negative least squares.
+
+    A fast non-negativity-constrained least squares
+    https://www.researchgate.net/publication/230554373_A_Fast_Non-negativity-constrained_Least_Squares_Algorithm
+    Rasmus Bro and Sijmen De Jong
+    Journal of Chemometrics. 11, 393–401 (1997)
+    """
+
+    n = len(A[0])
+
+    if not max_iters:
+        max_iters = n * 30
+
+    ATA = dot(transpose(A), A, dims=D2)
+    ATb = dot(transpose(A), b, dims=D2_D1)
+
+    x = zeros(n)  # type: Vector # type: ignore[assignment]
+    s = zeros(n)  # type: Vector # type: ignore[assignment]
+    w = subtract(ATb, dot(ATA, x, dims=D2_D1), dims=D1)  # type: Vector
+
+    # P tracks positive elements in x
+    P = [False] * n  # type: VectorBool
+
+    # Continue until all values of x are positive (non-negative results only)
+    # or we exhaust the iterations.
+    count = 0
+    while sum(P) < n and max(w[_i] for _i in range(n) if not P[_i]) > epsilon and count < max_iters:
+        # Find the index that maximizes w
+        # This will be an index not in P
+        imx = 0
+        mx = float('-inf')
+        for _i in range(n):
+            if not P[_i]:
+                temp = w[_i]
+                if temp > mx:
+                    imx = _i
+                    mx = temp
+        P[imx] = True
+
+        # Solve least squares problem for columns and rows not in P
+        idx = [_i for _i in range(n) if P[_i]]
+        v = dot(inv([[ATA[_i][_j] for _j in idx] for _i in idx]), [ATb[_i] for _i in idx], dims=D2_D1)
+        for _i, _v in zip(idx, v):
+            s[_i] = _v
+
+        # Deal with negative values
+        while _any([s[_i] <= epsilon for _i in range(n) if P[_i]]):
+            count += 1
+
+            # Calculate step size, alpha, to prevent any x from going negative
+            alpha = min(
+                [zdiv(x[_i], (x[_i] - s[_i]), float('inf')) for _i in range(n) if P[_i] * s[_i] <= epsilon]
+            )
+
+            # Update the solution
+            x = add(x, dot(alpha, subtract(s, x, dims=D1), dims=SC_D1), dims=D1)
+
+            # Remove indexes in P where x == 0
+            for _i in range(n):
+                if x[_i] <= epsilon:
+                    P[_i] = False
+
+            # Solve least squares problem again
+            idx = [_i for _i in range(n) if P[_i]]
+            v = dot(inv([[ATA[_i][_j] for _j in idx] for _i in idx]), [ATb[_i] for _i in idx], dims=D2_D1)
+            s = [0.0] * len(s)
+            for _i, _v in zip(idx, v):
+                s[_i] = _v
+
+        # Update the solution
+        x = s[:]
+        w = subtract(ATb, dot(ATA, x, dims=D2_D1), dims=D1)
+
+    # Return our final result, for better or for worse
+    res = math.hypot(*subtract(b, dot(A, x, dims=D2_D1), dims=D1))  # ||b-Ax||
+    return x, res

docs/src/dictionary/en-custom.txt

@@ -37,6 +37,7 @@ CVD
 CVDs
 Catmull
 Changelog
+Chemometrics
 Chroma
 Chromaticities
 Chromaticity
@@ -51,6 +52,7 @@ Cubehelix
 Culori
 Cz
 DCI
+De
 Deprecations
 Deregister
 Deregistration
@@ -95,6 +97,7 @@ IPT
 ITP
 ITU
 IgPgTg
+Illuminant
 Interpolator
 Itten
 Iz
@@ -105,6 +108,7 @@ JMh
 JND
 JSON
 Jacobian
+Jong
 Jsh
 Judd
 Jupyter
@@ -113,10 +117,13 @@ JzCzhz
 JzMzhz
 Jzazbz
 Kries
+Kubelka
 Kz
 LCh
 LChish
 LChuv
+LHTSS
+LLSS
 LMS
 Lab
 Labish
@@ -135,6 +142,7 @@ MkDocs
 Mollon
 Monochromacy
 Moroney
+Munk
 Mz
 NONINFRINGEMENT
 NaN
@@ -153,6 +161,7 @@ Oklrab
 Ostrowski
 Ottosson
 PQ
+Penrose
 Perceptibility
 Piecewise
 Planckian
@@ -177,6 +186,7 @@ RLAB
 ROMM
 RYB
 Raphson
+Rasmus
 SCD
 SDR
 SL
@@ -186,6 +196,7 @@ SVG
 Safdar
 Scalable
 Sharma
+Sijmen
 Sz
 TODO
 TORTIOUS
@@ -210,6 +221,7 @@ Vos
 Vz
 WCAG
 WCG
+Wijnen
 Wz
 XD
 XYB
@@ -256,16 +268,19 @@ desaturated
 deuteranomaly
 deuteranopia
 dichromacy
+differencing
 diffuser
 discretized
 docstring
 dyad
 easings
+emissive
 fixup
 formatter
 grayscale
 helixes
 hz
+illum
 illuminance
 illuminant
 illuminants
@@ -290,6 +305,9 @@ monochromacy
 monotonicity
 nd
 nm
+normalizations
+nx
+nxn
 oRGB
 opRGB
 opto
@@ -316,6 +334,7 @@ quantized
 quantizer
 rc
 reflectance
+reflectances
 repurpose
 rgb
 sRGB

pyproject.toml

@@ -105,23 +105,24 @@ lint.select = [
 ]
 
 lint.ignore = [
-    "E741",
+    "B905",
     "D202",
-    "D401",
-    "D212",
     "D203",
-    "N802",
+    "D212",
+    "D401",
+    "E741",
     "N801",
+    "N802",
     "N803",
     "N806",
     "N818",
-    "RUF012",
-    "RUF005",
     "PGH004",
+    "RUF002",
+    "RUF005",
+    "RUF012",
     "RUF022",
     "RUF023",
-    "RUF100",
-    "B905"
+    "RUF100"
 ]
 
 [tool.coverage.report]

tests/test_algebra.py

@@ -2613,6 +2613,59 @@ def dx2(x):
         # Ostrowski
         self.assertEqual(alg.solve_newton(1, f0, dx, ostrowski=True), (0.5, True))
 
+    def test_pinv(self):
+        """Test Moore-Penrose pseudo inverse."""
+
+        m = [
+            [0.4123907992659593, 0.3575843393838777, 0.1804807884018343],
+            [0.21263900587151033, 0.7151686787677553, 0.07219231536073373],
+            [0.019330818715591832, 0.11919477979462595, 0.9505321522496605]
+        ]
+
+        v = [0.047770200571454854, 0.02780940276126581, 0.22476064520055364]
+
+        # Negative results can be returned
+        result = alg.dot(alg.pinv(m), v)
+        self.assertEqual(result, [-5.551115123125783e-17, 0.015208514422912689, 0.23455058216100527])
+
+        wide = alg.pinv([[4, 5], [3, 3], [9, 7]])
+        self.assertEqual(alg.dot(wide, [3, 5, 6]), [0.29640718562873936, 0.538922155688625])
+
+        tall = alg.pinv([[4, 5, 3], [9, 7, 3]])
+        self.assertEqual(alg.dot(tall, [3, 5]), [0.2872727272727278, 0.2818181818181821, 0.14727272727272767])
+
+        with self.assertRaises(ValueError):
+            alg.pinv([1, 2, 3])
+
+    def test_fnnls(self):
+        """Test fast non-negative least squares method."""
+
+        m = [
+            [0.4123907992659593, 0.3575843393838777, 0.1804807884018343],
+            [0.21263900587151033, 0.7151686787677553, 0.07219231536073373],
+            [0.019330818715591832, 0.11919477979462595, 0.9505321522496605]
+        ]
+
+        v = [0.047770200571454854, 0.02780940276126581, 0.22476064520055364]
+
+        res = alg.fnnls(m, v)
+        b = alg.dot(alg.pinv(m), v)
+
+        # We should have no negative values, but we should be close to the `pinv` approach.
+        self.assertTrue(all(_a >= 0 for _a in res[0]))
+        self.assertTrue(res[1] < 1e-10)
+        self.assertTrue(all(math.isclose(_a, _b, rel_tol=1e-10, abs_tol=1e-11) for _a, _b in zip(res[0], b)))
+
+        # This is purposely beyond the range of a reasonable solution
+        # There will be residual
+        v = [0.6369580483012911, 0.262700212011267, 4.994106574466074e-17]
+        res = alg.fnnls(m, v)
+
+        # We should have no negative values, but we will have residual
+        self.assertFalse(res[1] < 1e-10)
+        self.assertTrue(all(_a >= 0 for _a in res[0]))
+        self.assertEqual(res[0], [1.477061311287275, 0.0, 0.0])
+
 
 def test_pprint(capsys):
     """Test matrix print."""