Repository: secretflow/spu (language: python)

Issue:
使用 SPU 实现非负矩阵分解算法(NMF)
> 此 ISSUE 为 [隐语开源共建计划（SecretFlow Open Source Contribution Plan，简称 SF OSCP）第二期](https://github.com/orgs/secretflow/discussions/741)任务 ISSUE，欢迎社区开发者参与共建～
> 若有感兴趣想要认领的任务，但还未报名，辛苦先完成[报名](https://www.wjx.top/vm/QKio0dq.aspx# )进行哈～


## 任务介绍

- 任务名称：使用 SPU 实现非负矩阵分解算法(NMF)
- 技术方向：SPU/SML
- 任务难度：挑战🌟🌟🌟

## 详细要求
- 安全性（尽量少 reveal）
- 功能性：实现算法的基本功能，包括：
  - 支持 fit, transform, inverse_transform
  - 同时支持 element-wise L1 norm 和 Frobenius norm
- 收敛性：包含 simulator 跑出的实验数据并且证明收敛性/准确性（最好有与明文 sklearn 结果的对比，可详见[参考 PR](https://github.com/secretflow/spu/pull/240)）
- 正确性：请确保提交的代码内容为可以直接运行的
- 代码规范：Python 代码需要使用 black+isort 进行格式化（流水线包含代码规范检查卡点）
- 提交说明：关联该 isuue 并提交代码至 https://github.com/secretflow/spu/tree/main/sml
- 特殊说明：若某个特性有特殊的限制，如需要 FM128，需要更多 fxp 等需要在注释文档中明确说明

## 能力要求
- 熟悉经典的机器学习算法
- 熟悉 JAX 或 NumPy，可以使用 NumPy 实现算法

## 操作说明
- 操作指引：https://www.secretflow.org.cn/docs/spu/latest/en-US/getting_started/tutorials/develop_your_first_mpc_application
- 可参考明文实现：[sklearn.decomposition](https://scikit-learn.org/stable/modules/generated/sklearn.decomposition.NMF.html#sklearn.decomposition.NMF)
-  [范例](https://github.com/secretflow/spu/blob/main/sml/linear_model/simple_sgd.py)
-  [参考 PR](https://github.com/secretflow/spu/pull/240)

## 开发须知
以下部分代码请必须增加代码注释，对对应代码模块进行说明，包括：
-  __init__函数的超参数含义
-  fit 的具体算法实现说明


Discussion / hints:
tarantula-leo give it to me.

Fix patch:
diff --git a/sml/decomposition/BUILD.bazel b/sml/decomposition/BUILD.bazel
index 82297199..c80a6751 100644
--- a/sml/decomposition/BUILD.bazel
+++ b/sml/decomposition/BUILD.bazel
@@ -21,3 +21,8 @@ py_library(
     srcs = ["pca.py"],
     deps = ["//sml/utils:extmath"],
 )
+
+py_library(
+    name = "nmf",
+    srcs = ["nmf.py"],
+)
diff --git a/sml/decomposition/emulations/BUILD.bazel b/sml/decomposition/emulations/BUILD.bazel
index 18807994..100d6e40 100644
--- a/sml/decomposition/emulations/BUILD.bazel
+++ b/sml/decomposition/emulations/BUILD.bazel
@@ -35,6 +35,16 @@ py_binary(
     ],
 )
 
+py_binary(
+    name = "nmf_emul",
+    srcs = ["nmf_emul.py"],
+    data = [":conf"],
+    deps = [
+        "//sml/decomposition:nmf",
+        "//sml/utils:emulation",
+    ],
+)
+
 filegroup(
     name = "conf",
     srcs = [
diff --git a/sml/decomposition/emulations/nmf_emul.py b/sml/decomposition/emulations/nmf_emul.py
new file mode 100644
index 00000000..3304312e
--- /dev/null
+++ b/sml/decomposition/emulations/nmf_emul.py
@@ -0,0 +1,154 @@
+# Copyright 2023 Ant Group Co., Ltd.
+#
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+#
+#      https://www.apache.org/licenses/LICENSE-2.0
+#
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+
+import os
+import sys
+
+import numpy as np
+from sklearn.decomposition import NMF as SklearnNMF
+
+# Add the sml directory to the path
+sys.path.append(os.path.join(os.path.dirname(__file__), '../../../'))
+
+import sml.utils.emulation as emulation
+from sml.decomposition.nmf import NMF
+
+
+def test_nmf(mode: emulation.Mode.MULTIPROCESS):
+    def emul_nmf():
+        def proc1(X, random_matrixA, random_matrixB):
+            model = NMF(
+                n_components=n_components,
+                l1_ratio=l1_ratio,
+                alpha_W=alpha_W,
+                random_matrixA=random_matrixA,
+                random_matrixB=random_matrixB,
+            )
+
+            W = model.fit_transform(X)
+            H = model._components
+            X_reconstructed = model.inverse_transform(W)
+            err = model.reconstruction_err_
+            return W, H, X_reconstructed, err
+
+        # Create a simple dataset and random_matrix
+        X = np.random.randint(1, 100, (1000, 10))
+        X = np.array(X, dtype=float)
+        n_samples, n_features = X.shape
+        n_components = 5
+        random_seed = 0
+        random_state = np.random.RandomState(random_seed)
+        A = random_state.standard_normal(size=(n_components, n_features))
+        B = random_state.standard_normal(size=(n_samples, n_components))
+        l1_ratio = 0.1
+        alpha_W = 0.01
+
+        X_spu = emulator.seal(X)
+        A_spu = emulator.seal(A)
+        B_spu = emulator.seal(B)
+
+        # Run the emulation
+        W, H, X_reconstructed, err = emulator.run(proc1)(X_spu, A_spu, B_spu)
+        print("reconstruction_error: ", err)
+
+        # sklearn
+        model = SklearnNMF(
+            n_components=n_components,
+            init='random',
+            random_state=random_seed,
+            l1_ratio=l1_ratio,
+            solver="mu",
+            alpha_W=alpha_W,
+        )
+        W_Sklearn = model.fit_transform(X)
+        H_Sklearn = model.components_
+        X_reconstructed_Sklearn = model.inverse_transform(W_Sklearn)
+        err = model.reconstruction_err_
+        print("reconstruction_error_sklearn: ", err)
+        assert np.allclose(W_Sklearn, W, atol=5e-1)
+        assert np.allclose(H_Sklearn, H, atol=5e-1)
+        assert np.allclose(X_reconstructed_Sklearn, X_reconstructed, atol=5e-1)
+
+    def emul_nmf_seperate():
+        def proc2(X, random_matrixA, random_matrixB):
+            model = NMF(
+                n_components=n_components,
+                l1_ratio=l1_ratio,
+                alpha_W=alpha_W,
+                random_matrixA=random_matrixA,
+                random_matrixB=random_matrixB,
+            )
+
+            model.fit(X)
+            W = model.transform(X, transform_iter=40)
+            H = model._components
+            X_reconstructed = model.inverse_transform(W)
+            return W, H, X_reconstructed
+
+        # Create a simple dataset and random_matrix
+        X = np.random.randint(1, 100, (1000, 10))
+        X = np.array(X, dtype=float)
+        n_samples, n_features = X.shape
+        n_components = 5
+        random_seed = 0
+        random_state = np.random.RandomState(random_seed)
+        A = random_state.standard_normal(size=(n_components, n_features))
+        B = random_state.standard_normal(size=(n_samples, n_components))
+        l1_ratio = 0.1
+        alpha_W = 0.01
+
+        X_spu = emulator.seal(X)
+        A_spu = emulator.seal(A)
+        B_spu = emulator.seal(B)
+
+        # Run the emulation_seperate
+        W_seperate, H_seperate, X_reconstructed_seperate = emulator.run(proc2)(
+            X_spu, A_spu, B_spu
+        )
+        print("W_matrix_spu_seperate: ", W_seperate[:5, :5])
+        print("H_matrix_spu_seperate: ", H_seperate[:5, :5])
+        print("X_reconstructed_spu_seperate: ", X_reconstructed_seperate[:5, :5])
+
+        # sklearn_seperate
+        model = SklearnNMF(
+            n_components=n_components,
+            init='random',
+            random_state=random_seed,
+            l1_ratio=l1_ratio,
+            solver="mu",
+            alpha_W=alpha_W,
+        )
+        model.fit(X)
+        W_Sklearn_seperate = model.transform(X)
+        H_Sklearn_seperate = model.components_
+        X_reconstructed_Sklearn_seperate = model.inverse_transform(W_Sklearn_seperate)
+        print("W_matrix_sklearn_seperate: ", W_Sklearn_seperate[:5, :5])
+        print("H_matrix_sklearn_seperate: ", H_Sklearn_seperate[:5, :5])
+        print(
+            "X_reconstructed_sklearn_seperate: ",
+            X_reconstructed_Sklearn_seperate[:5, :5],
+        )
+
+    try:
+        conf_path = "sml/decomposition/emulations/3pc.json"
+        emulator = emulation.Emulator(conf_path, mode, bandwidth=300, latency=20)
+        emulator.up()
+        emul_nmf()
+        emul_nmf_seperate()
+    finally:
+        emulator.down()
+
+
+if __name__ == "__main__":
+    test_nmf(emulation.Mode.MULTIPROCESS)
diff --git a/sml/decomposition/nmf.py b/sml/decomposition/nmf.py
new file mode 100644
index 00000000..680c3483
--- /dev/null
+++ b/sml/decomposition/nmf.py
@@ -0,0 +1,222 @@
+# Copyright 2023 Ant Group Co., Ltd.
+#
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+#
+#      https://www.apache.org/licenses/LICENSE-2.0
+#
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+
+import numpy as np
+import jax.numpy as jnp
+
+
+def update_w(X, W, H, H_sum, HHt, XHt, l1_reg_W, l2_reg_W):
+    # H_sum is not used in 'frobenius'
+    if XHt is None:
+        XHt = jnp.dot(X, H.T)
+    numerator = XHt
+    if HHt is None:
+        HHt = jnp.dot(H, H.T)
+    denominator = jnp.dot(W, HHt)
+    denominator = denominator + l1_reg_W + l2_reg_W * W
+    numerator /= denominator
+    delta_W = numerator
+    W *= delta_W
+    return W, H_sum, HHt, XHt
+
+
+def update_h(X, W, H, l1_reg_H, l2_reg_H):
+    numerator = jnp.dot(W.T, X)
+    denominator = jnp.linalg.multi_dot([W.T, W, H])
+    denominator = denominator + l1_reg_H + l2_reg_H * H
+    delta_H = numerator
+    delta_H /= denominator
+    H *= delta_H
+    return H
+
+
+def init_nmf(X, n_components, random_matrixA, random_matrixB):
+    avg = jnp.sqrt(np.mean(X) / n_components)
+    W = avg * random_matrixB.astype(X.dtype)
+    W_new = jnp.abs(W)
+    H = avg * random_matrixA.astype(X.dtype)
+    H_new = jnp.abs(H)
+    return W_new, H_new
+
+
+def compute_regularization(X, l1_ratio, alpha_W, alpha_H):
+    n_samples, n_features = X.shape
+    alpha_H = alpha_W
+    l1_reg_W = n_features * alpha_W * l1_ratio
+    l1_reg_H = n_samples * alpha_H * l1_ratio
+    l2_reg_W = n_features * alpha_W * (1.0 - l1_ratio)
+    l2_reg_H = n_samples * alpha_H * (1.0 - l1_ratio)
+    return l1_reg_W, l1_reg_H, l2_reg_W, l2_reg_H
+
+
+def _beta_divergence(X, W, H, beta, square_root=False):
+    if beta == 'frobenius':
+        vec = jnp.ravel(X - jnp.dot(W, H))
+        res = jnp.dot(vec, vec) / 2.0
+        if square_root:
+            res = jnp.sqrt(res * 2)
+    return res
+
+
+class NMF:
+    def __init__(
+        self,
+        n_components: int,
+        max_iter: int = 200,
+        l1_ratio: float = 0.0,
+        alpha_W: float = 0.0,
+        beta_loss: str = 'frobenius',
+        alpha_H=None,
+        random_matrixA=None,
+        random_matrixB=None,
+    ):
+        """Compute Non-negative Matrix Factorization (NMF).
+
+        Find two non-negative matrices (W, H) whose product approximates the non-
+        negative matrix X.
+
+        Parameters
+        ----------
+        n_components : int, default=None
+            Number of components.
+        max_iter : int, default=200
+            Maximum number of iterations.
+        l1_ratio : float, default=0.0
+            The regularization mixing parameter, with 0 <= l1_ratio <= 1.
+            For l1_ratio = 0 the penalty is an elementwise L2 penalty
+            (aka Frobenius Norm).
+            For l1_ratio = 1 it is an elementwise L1 penalty.
+            For 0 < l1_ratio < 1, the penalty is a combination of L1 and L2.
+        alpha_W : float, default=0.0
+            Constant that multiplies the regularization terms of `W`. Set it to zero
+            (default) to have no regularization on `W`.
+        alpha_H : None
+            Constant that multiplies the regularization terms of `H`. Set it takes
+            the same value as `alpha_W`.
+        beta_loss : str, default='frobenius'
+            Beta divergence to be minimized, measuring the distance between X
+            and the dot product WH. Only support 'frobenius' now.
+        random_matrixA : None
+            Non-negative random matrices, scaled with:
+            sqrt(X.mean() / n_components), represents matrix H
+        random_matrixB : None
+            Non-negative random matrices, scaled with:
+            sqrt(X.mean() / n_components), represents matrix W
+
+        References
+        ----------
+        Lee, D. D., & Seung, H., S. (2001). Algorithms for Non-negative Matrix
+        Factorization. Adv. Neural Inform. Process. Syst.. 13.
+        """
+        # parameter check.
+        assert n_components > 0, f"n_components should >0"
+        assert beta_loss == 'frobenius', f"beta_loss only support frobenius now"
+
+        if alpha_H is None:
+            alpha_H = alpha_W
+        self._beta_loss = beta_loss
+        self._n_components = n_components
+        self._max_iter = max_iter
+        self._l1_ratio = l1_ratio
+        self._alpha_W = alpha_W
+        self._alpha_H = alpha_H
+        self._random_matrixA = random_matrixA
+        self._random_matrixB = random_matrixB
+        self._update_H = True
+        self._components = None
+
+    def fit(self, X):
+        """Learn a NMF model for the data X.
+        In the 'mu' solver, we use the multiplicative update algorithm to compute matrix W and H.
+        Parameters
+        ----------
+        X : {array-like}, shape (n_samples, n_features)
+            Training data.
+        Returns
+        -------
+        self : object
+            Returns an instance of self.
+        Notes
+        -----
+        (1) To prevent overflow error when using large data sets or get more accurate results,
+        you can modify the definition of simulator as follows:
+        config = spu_pb2.RuntimeConfig(
+            protocol=spu_pb2.ProtocolKind.ABY3,
+            field=spu_pb2.FieldType.FM128,
+            fxp_fraction_bits=30,
+            )
+        sim_aby = spsim.Simulator(3, config)
+        (2) The regularization parameter is deprecated in sklearn_1.0 and is removed in sklearn_1.2.
+        Use 'alpha_W' and 'alpha_H' instead.
+        (3) It's not very convenient to use jax_control_flow in spu right now, so we use 'max_iter' to control
+        the iteration stop, and you can use 'reconstruction_err_' to learn about Frobenius norm of the matrix difference.
+        If you want to get closer to the sklearn, need to keep the number of iterations consistent.
+        (4) Only support 'random' init and 'mu' solver now.
+        """
+        self.fit_transform(X)
+        return self
+
+    def transform(self, X, transform_iter=40):
+        assert self._components is not None, f"should fit before transform"
+        self._update_H = False
+        self._max_iter = transform_iter
+        W = self.fit_transform(X)
+        return W
+
+    def fit_transform(self, X):
+        # check random_matrix shape
+        assert self._random_matrixA.shape == (
+            self._n_components,
+            X.shape[1],
+        ), f"Expected random_matrixA to be ({self._n_components}, {X.shape[1]}) array, got {self._random_matrixA.shape}"
+        assert self._random_matrixB.shape == (
+            X.shape[0],
+            self._n_components,
+        ), f"Expected random_matrixA to be ({X.shape[0]}, {self._n_components}) array, got {self._random_matrixB.shape}"
+
+        # init matrix W&H
+        if self._update_H:
+            W, H = init_nmf(
+                X, self._n_components, self._random_matrixA, self._random_matrixB
+            )
+        else:
+            avg = jnp.sqrt(X.mean() / self._n_components)
+            W = jnp.full((X.shape[0], self._n_components), avg, dtype=X.dtype)
+            H = self._components
+
+        # compute the regularization parameters
+        l1_reg_W, l1_reg_H, l2_reg_W, l2_reg_H = compute_regularization(
+            X, self._l1_ratio, self._alpha_W, self._alpha_H
+        )
+        H_sum, HHt, XHt = None, None, None
+
+        # use multiplicative update solver
+        for _ in range(self._max_iter):
+            W, H_sum, HHt, XHt = update_w(X, W, H, H_sum, HHt, XHt, l1_reg_W, l2_reg_W)
+            if self._update_H:
+                H = update_h(X, W, H, l1_reg_H, l2_reg_H)
+                H_sum, HHt, XHt = None, None, None
+
+        # compute the reconstruction error
+        if self._update_H:
+            self._components = H
+            self.reconstruction_err_ = _beta_divergence(
+                X, W, H, self._beta_loss, square_root=True
+            )
+
+        return W
+
+    def inverse_transform(self, X):
+        assert self._components is not None, f"should fit before inverse_transform"
+        return jnp.dot(X, self._components)


Test patch:
diff --git a/sml/decomposition/tests/BUILD.bazel b/sml/decomposition/tests/BUILD.bazel
index b785cbe8..fff72220 100644
--- a/sml/decomposition/tests/BUILD.bazel
+++ b/sml/decomposition/tests/BUILD.bazel
@@ -25,3 +25,13 @@ py_test(
         "//spu/utils:simulation",
     ],
 )
+
+py_test(
+    name = "nmf_test",
+    srcs = ["nmf_test.py"],
+    deps = [
+        "//sml/decomposition:nmf",
+        "//spu:init",
+        "//spu/utils:simulation",
+    ],
+)
diff --git a/sml/decomposition/tests/nmf_test.py b/sml/decomposition/tests/nmf_test.py
new file mode 100644
index 00000000..44ec7b8e
--- /dev/null
+++ b/sml/decomposition/tests/nmf_test.py
@@ -0,0 +1,157 @@
+# Copyright 2023 Ant Group Co., Ltd.
+#
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+#
+#      https://www.apache.org/licenses/LICENSE-2.0
+#
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+
+import os
+import sys
+import unittest
+
+import numpy as np
+from sklearn.decomposition import NMF as SklearnNMF
+import spu.spu_pb2 as spu_pb2
+import spu.utils.simulation as spsim
+
+# Add the sml directory to the path
+sys.path.append(os.path.join(os.path.dirname(__file__), '../../../'))
+
+from sml.decomposition.nmf import NMF
+
+
+class UnitTests(unittest.TestCase):
+    def test_nmf(self):
+        config = spu_pb2.RuntimeConfig(
+            protocol=spu_pb2.ProtocolKind.ABY3,
+            field=spu_pb2.FieldType.FM128,
+            fxp_fraction_bits=30,
+        )
+        sim = spsim.Simulator(3, config)
+
+        # Test fit_transform
+        def proc1(X, random_matrixA, random_matrixB):
+            model = NMF(
+                n_components=n_components,
+                l1_ratio=l1_ratio,
+                alpha_W=alpha_W,
+                random_matrixA=random_matrixA,
+                random_matrixB=random_matrixB,
+            )
+
+            W = model.fit_transform(X)
+            H = model._components
+            X_reconstructed = model.inverse_transform(W)
+            err = model.reconstruction_err_
+            return W, H, X_reconstructed, err
+
+        # Create a simple dataset and random_matrix
+        X = np.random.randint(1, 100, (1000, 10))
+        X = np.array(X, dtype=float)
+        n_samples, n_features = X.shape
+        n_components = 5
+        random_seed = 0
+        random_state = np.random.RandomState(random_seed)
+        A = random_state.standard_normal(size=(n_components, n_features))
+        B = random_state.standard_normal(size=(n_samples, n_components))
+        l1_ratio = 0.1
+        alpha_W = 0.01
+
+        # Run the simulation
+        W, H, X_reconstructed, err = spsim.sim_jax(sim, proc1)(X, A, B)
+        print("reconstruction_error: ", err)
+
+        # sklearn
+        model = SklearnNMF(
+            n_components=n_components,
+            init='random',
+            random_state=random_seed,
+            l1_ratio=l1_ratio,
+            solver="mu",
+            alpha_W=alpha_W,
+        )
+        W_Sklearn = model.fit_transform(X)
+        H_Sklearn = model.components_
+        X_reconstructed_Sklearn = model.inverse_transform(W_Sklearn)
+        err = model.reconstruction_err_
+        print("reconstruction_error_sklearn: ", err)
+        self.assertTrue(np.allclose(W_Sklearn, W, atol=5e-1))
+        self.assertTrue(np.allclose(H_Sklearn, H, atol=5e-1))
+        self.assertTrue(
+            np.allclose(X_reconstructed_Sklearn, X_reconstructed, atol=5e-1)
+        )
+
+    def test_nmf_seperate(self):
+        config = spu_pb2.RuntimeConfig(
+            protocol=spu_pb2.ProtocolKind.ABY3,
+            field=spu_pb2.FieldType.FM128,
+            fxp_fraction_bits=30,
+        )
+        sim = spsim.Simulator(3, config)
+
+        # Test fit and transform
+        def proc2(X, random_matrixA, random_matrixB):
+            model = NMF(
+                n_components=n_components,
+                l1_ratio=l1_ratio,
+                alpha_W=alpha_W,
+                random_matrixA=random_matrixA,
+                random_matrixB=random_matrixB,
+            )
+
+            model.fit(X)
+            W = model.transform(X, transform_iter=40)
+            H = model._components
+            X_reconstructed = model.inverse_transform(W)
+            return W, H, X_reconstructed
+
+        # Create a simple dataset and random_matrix
+        X = np.random.randint(1, 100, (1000, 10))
+        X = np.array(X, dtype=float)
+        n_samples, n_features = X.shape
+        n_components = 5
+        random_seed = 0
+        random_state = np.random.RandomState(random_seed)
+        A = random_state.standard_normal(size=(n_components, n_features))
+        B = random_state.standard_normal(size=(n_samples, n_components))
+        l1_ratio = 0.1
+        alpha_W = 0.01
+
+        # Run the simulation_seperate
+        W_seperate, H_seperate, X_reconstructed_seperate = spsim.sim_jax(sim, proc2)(
+            X, A, B
+        )
+        print("W_matrix_spu_seperate: ", W_seperate[:5, :5])
+        print("H_matrix_spu_seperate: ", H_seperate[:5, :5])
+        print("X_reconstructed_spu_seperate: ", X_reconstructed_seperate[:5, :5])
+
+        # sklearn_seperate
+        model = SklearnNMF(
+            n_components=n_components,
+            init='random',
+            random_state=random_seed,
+            l1_ratio=l1_ratio,
+            solver="mu",
+            alpha_W=alpha_W,
+        )
+        model.fit(X)
+        W_Sklearn_seperate = model.transform(X)
+        H_Sklearn_seperate = model.components_
+        X_reconstructed_Sklearn_seperate = model.inverse_transform(W_Sklearn_seperate)
+        print("W_matrix_sklearn_seperate: ", W_Sklearn_seperate[:5, :5])
+        print("H_matrix_sklearn_seperate: ", H_Sklearn_seperate[:5, :5])
+        print(
+            "X_reconstructed_sklearn_seperate: ",
+            X_reconstructed_Sklearn_seperate[:5, :5],
+        )
+
+
+if __name__ == "__main__":
+    unittest.main()