rpca pcp now provides LQ decomposition

qolmat/imputations/rpca/rpca_pcp.py

@@ -80,14 +80,18 @@ def get_params_scale(self, D: NDArray):
         dict_params = {"mu": mu, "lam": lam}
         return dict_params
 
-    def decompose_on_basis(self, D: NDArray, Omega: NDArray, Q: NDArray) -> NDArray:
+    def decompose_on_basis(
+        self, D: NDArray, Omega: NDArray, Q: NDArray
+    ) -> Tuple[NDArray, NDArray]:
         n_rows, n_cols = D.shape
         if n_rows == 1 or n_cols == 1:
             return D, np.full_like(D, 0)
         M, A, L, Q = self.decompose_rpca(D, Omega)
         return M, A
 
-    def decompose_rpca(self, D: NDArray, Omega: NDArray) -> Tuple[NDArray, NDArray, None, None]:
+    def decompose_rpca(
+        self, D: NDArray, Omega: NDArray
+    ) -> Tuple[NDArray, NDArray, NDArray, NDArray]:
         """
         Estimate the relevant parameters then compute the PCP RPCA decomposition
 
@@ -123,7 +127,8 @@ def decompose_rpca(self, D: NDArray, Omega: NDArray) -> Tuple[NDArray, NDArray,
 
         M: NDArray = D - A
         for iteration in range(self.max_iterations):
-            M = rpca_utils.svd_thresholding(D - A + Y / mu, 1 / mu)
+            L, Q = rpca_utils.svd_thresholding(D - A + Y / mu, 1 / mu)
+            M = L @ Q
             A = rpca_utils.soft_thresholding(D - M + Y / mu, lam / mu)
             A[~Omega] = (D - M)[~Omega]
 
@@ -137,7 +142,7 @@ def decompose_rpca(self, D: NDArray, Omega: NDArray) -> Tuple[NDArray, NDArray,
 
         self._check_cost_function_minimized(D, M, A, Omega, lam)
 
-        return M, A, None, None
+        return M, A, L, Q
 
     def _check_cost_function_minimized(
         self,

qolmat/imputations/rpca/rpca_utils.py

@@ -2,6 +2,7 @@
 Modular utility functions for RPCA
 """
 
+from typing import Tuple
 import numpy as np
 from numpy.typing import NDArray
 import scipy
@@ -60,7 +61,7 @@ def soft_thresholding(
     return np.sign(X) * np.maximum(np.abs(X) - threshold, 0)
 
 
-def svd_thresholding(X: NDArray, threshold: float) -> NDArray:
+def svd_thresholding(X: NDArray, threshold: float) -> Tuple[NDArray, NDArray]:
     """
     Apply the shrinkage operator to the singular values obtained from the SVD of X.
 
@@ -73,16 +74,20 @@ def svd_thresholding(X: NDArray, threshold: float) -> NDArray:
 
     Returns
     -------
-    NDArray
-        Array obtained by computing U * shrink(s) * V where
+    Tuple[NDArray, NDArray]
+        Two arrays L and Q of minimal Frobenius norm such that L @ Q = U * shrink(s) * V where
             U is the array of left singular vectors of X
             V is the array of the right singular vectors of X
             s is the array of the singular values as a diagonal array
+        L and Q minimize
     """
 
     U, s, Vh = np.linalg.svd(X, full_matrices=False)
     s = soft_thresholding(s, threshold)
-    return U @ (np.diag(s) @ Vh)
+    # return U @ (np.diag(s) @ Vh)
+    L = U @ np.sqrt(np.diag(s))
+    Q = np.sqrt(np.diag(s)) @ Vh
+    return L, Q
 
 
 def l1_norm(M: NDArray) -> float:

tests/imputations/rpca/test_rpca_utils.py

@@ -55,8 +55,8 @@ def test_rpca_utils_soft_thresholding(X: NDArray, threshold: float):
 @pytest.mark.parametrize("X", [X_complete])
 @pytest.mark.parametrize("threshold", [0.95])
 def test_rpca_utils_svd_thresholding(X: NDArray, threshold: float):
-    result = svd_thresholding(X=X, threshold=threshold)
-    print(result)
+    L_result, Q_result = svd_thresholding(X=X, threshold=threshold)
+    result = L_result @ Q_result
     X_expected = np.array(
         [
             [0.928, 6.182, 3.857, 3.857],
@@ -81,6 +81,12 @@ def test_rpca_utils_toeplitz_matrix(T: int, dimension: int):
     result = toeplitz_matrix(T=T, dimension=dimension)
     result_np = result.toarray()
     X_exp = np.array(
-        [[1, 0, -1, 0, 0], [0, 1, 0, -1, 0], [0, 0, 1, 0, -1], [0, 0, 0, 0, 0], [0, 0, 0, 0, 0]]
+        [
+            [1, 0, -1, 0, 0],
+            [0, 1, 0, -1, 0],
+            [0, 0, 1, 0, -1],
+            [0, 0, 0, 0, 0],
+            [0, 0, 0, 0, 0],
+        ]
     )
     np.testing.assert_allclose(result_np, X_exp)