Repository: secretflow/spu (language: python)

Issue:
[Perceptron-sgd]使用 SPU 实现感知机(perceptron)算法
> 此 ISSUE 为 [隐语开源共建计划（SecretFlow Open Source Contribution Plan，简称 SF OSCP）第二期](https://github.com/orgs/secretflow/discussions/741) [使用 SPU 实现感知机(perceptron)算法](https://github.com/secretflow/spu/issues/250) 子任务，若对该任务感兴趣，欢迎在 https://github.com/secretflow/spu/issues/250 提交你的不同设计思路。
> 也欢迎其他开发者点击[这里](https://github.com/orgs/secretflow/discussions/741)，选择其他你感兴趣的任务，热烈欢迎社区开发者参与共建～

## 设计思路
by：[zoeStartover](https://github.com/zoeStartover)
1. 初始化感知机模型，确定是否需要正则(L1,L2,ElasticNet,None)，是否需要拟合bias项等基本参数（参照sklearn）；输入训练集X,Y = (x_i,y_i)，优化目的是最小化-sum(y_i(w*x_i+b))，每一圈使用随机梯度下降法(SGD)更新参数w,b，w=w+η * y_i * x_i，b=b+ η * y_i，其中可选择的正则化方法包括L1,L2,ElasticNet(L1+L2)。
2. fit方法：训练感知机模型，输入训练集对X,Y = (x_i,y_i)，若某点计算为误分类点，更新参数w和b，迭代上述过程至max_iter；
3. predict方法：输入待预测样本x，输出预测结果y=sign(w*x+b)。


## 任务介绍

- 任务名称：使用 SPU 实现感知机(perceptron)算法
- 技术方向：SPU/SML
- 任务难度：热身🌟

## 详细要求
- 安全性（尽量少 reveal）
- 功能性：实现算法的基本功能，包括：
  - 支持 fit，predict
  - 提供配置是否需要拟合 bias 项，是否需要正则等基本参数
  - 支持起码一种正则方式（L2, L1, elastic...）
- 收敛性：包含 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
- 可参考明文实现：[linear_model/_perceptron](https://github.com/scikit-learn/scikit-learn/blob/7f9bad99d/sklearn/linear_model/_perceptron.py#L9)
-  [范例](https://github.com/secretflow/spu/blob/main/sml/linear_model/simple_sgd.py)
-  [参考 PR](https://github.com/secretflow/spu/pull/240)


## 认领说明
本任务可有多种实现方式，故支持一个任务有多位开发者进行认领，请在认领任务后，在该 issue 下 comment 你的具体设计思路。

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


Discussion / hints:
zoeStartover Give it to me
1. 已使用jax实现perceptron的fit,predict（包含正则、是否需要拟合bias等多个参数，SGD优化）；
2. 使用鸢尾花数据集（100 samples）与sklearn对比准确率差1%；
3. 已使用 black+isort 进行格式化。

Fix patch:
diff --git a/sml/linear_model/BUILD.bazel b/sml/linear_model/BUILD.bazel
index 9c3e564e..53b52d6c 100644
--- a/sml/linear_model/BUILD.bazel
+++ b/sml/linear_model/BUILD.bazel
@@ -36,3 +36,8 @@ py_library(
     name = "ridge",
     srcs = ["ridge.py"],
 )
+
+py_library(
+    name = "pla",
+    srcs = ["pla.py"],
+)
diff --git a/sml/linear_model/emulations/BUILD.bazel b/sml/linear_model/emulations/BUILD.bazel
index 8a6c5245..b691eb56 100644
--- a/sml/linear_model/emulations/BUILD.bazel
+++ b/sml/linear_model/emulations/BUILD.bazel
@@ -44,3 +44,12 @@ py_binary(
         "//sml/utils:emulation",
     ],
 )
+
+py_binary(
+    name = "pla_emul",
+    srcs = ["pla_emul.py"],
+    deps = [
+        "//sml/linear_model:pla",
+        "//sml/utils:emulation",
+    ],
+)
diff --git a/sml/linear_model/emulations/pla_emul.py b/sml/linear_model/emulations/pla_emul.py
new file mode 100644
index 00000000..a602bcc4
--- /dev/null
+++ b/sml/linear_model/emulations/pla_emul.py
@@ -0,0 +1,102 @@
+# 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
+#
+#   http://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 jax.numpy as jnp
+import pandas as pd
+import sklearn.linear_model as sk
+from sklearn.datasets import load_iris
+
+import sml.utils.emulation as emulation
+from sml.linear_model.pla import Perceptron
+
+
+def emul_perceptron(mode: emulation.Mode.MULTIPROCESS):
+    def proc(x, y):
+        model = Perceptron(
+            max_iter=20,
+            eta0=1.0,
+            penalty='elasticnet',
+            alpha=0.001,
+            fit_intercept=True,
+            l1_ratio=0.7,
+            patience=10,
+            batch_size=64,
+            early_stop=True,
+        )
+
+        return model.fit(x, y).predict(x)
+
+    def load_data():
+        iris = load_iris()
+        df = pd.DataFrame(
+            iris.data,
+            columns=['sepal length', 'sepal width', 'petal length', 'petal width'],
+        )
+        df['label'] = iris.target
+
+        # only use sepal length and sepal width features
+        # 100 samples
+        data = jnp.array(df.iloc[0:100, [0, 1, -1]])
+        x, y = data[:, :-1], data[:, -1]
+
+        # y is -1 or 1
+        y = jnp.sign(y)
+        y = jnp.where(y <= 0, -1, y)
+        y = y.reshape((y.shape[0], 1))
+
+        return x, y
+
+    try:
+        # bandwidth and latency only work for docker mode
+        emulator = emulation.Emulator(
+            emulation.CLUSTER_ABY3_3PC, mode, bandwidth=300, latency=20
+        )
+        emulator.up()
+
+        # load mock data
+        x, y = load_data()
+        n_samples = len(y)
+
+        # compare with sklearn
+        sk_pla = sk.Perceptron(
+            max_iter=20,
+            eta0=1.0,
+            penalty='elasticnet',
+            alpha=0.001,
+            l1_ratio=0.7,
+            fit_intercept=True,
+        )
+        result_sk = sk_pla.fit(x, y).predict(x)
+        result_sk = result_sk.reshape(result_sk.shape[0], 1)
+        acc_sk = jnp.sum((result_sk == y)) / n_samples * 100
+
+        # mark these data to be protected in SPU
+        x_spu, y_spu = emulator.seal(x, y)
+
+        # run
+        result = emulator.run(proc)(x_spu, y_spu)
+        result = result.reshape(result.shape[0], 1)
+        acc_ = jnp.sum((result == y)) / n_samples * 100
+
+        # print acc
+        print(f"Accuracy in SKlearn: {acc_sk:.2f}%")
+        print(f"Accuracy in SPU: {acc_:.2f}%")
+
+    finally:
+        emulator.down()
+
+
+if __name__ == "__main__":
+    emul_perceptron(emulation.Mode.MULTIPROCESS)
diff --git a/sml/linear_model/pla.py b/sml/linear_model/pla.py
new file mode 100644
index 00000000..4acd450c
--- /dev/null
+++ b/sml/linear_model/pla.py
@@ -0,0 +1,242 @@
+# 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
+#
+#   http://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 math
+from enum import Enum
+
+import jax
+import jax.numpy as jnp
+
+
+class Penalty(Enum):
+    NONE = None
+    L1 = 'l1'
+    L2 = 'l2'
+    EN = 'elasticnet'
+
+
+class Perceptron:
+    """Perceptron.
+
+    Parameters
+    ----------
+
+    penalty : {'l2','l1','elasticnet',None(NoneType)}, default=None
+        The penalty (aka regularization term) to be used.
+
+    alpha : float, default=0.0001
+        Constant that multiplies the regularization term if regularization is
+        used.
+
+    patience: int, default=10
+        How long to wait after last time loss improved.
+
+    l1_ratio : float, default=0.15
+        The Elastic Net mixing parameter, with `0 <= l1_ratio <= 1`.
+        `l1_ratio=0` corresponds to L2 penalty, `l1_ratio=1` to L1.
+        Only used if `penalty='elasticnet'`.
+
+    fit_intercept : bool, default=True
+        Whether the intercept should be estimated or not. If False, the
+        data is assumed to be already centered.
+
+    max_iter : int, default=1000
+        The maximum number of passes over the training data (aka epochs).
+        It only impacts the behavior in the ``fit`` method, and not the
+        :meth:`partial_fit` method.
+
+    eta0 : float, default=1
+        Constant by which the updates are multiplied.
+
+    batch_size : int, default=-1
+        The batch size is a number of samples processed before the model is updated.
+        if batch_size = -1, all samples will be processed in an iteration.
+
+    early_stop : bool, default=True
+        When the loss function of the model no longer decreases, the early stop method
+        can avoid the overtraining of the model if True.
+    """
+
+    def __init__(
+        self,
+        penalty=None,
+        alpha=0.0001,
+        patience=10,
+        l1_ratio=0.15,
+        fit_intercept=True,
+        max_iter=1000,
+        eta0=1.0,
+        batch_size=-1,
+        early_stop=True,
+    ):
+        # parameter check.
+        assert max_iter > 0, f"max_iter should >0"
+        assert eta0 > 0, f"eta0 should >0"
+        assert alpha > 0, f"alpha should >0"
+        assert patience > 0, f"patience should >0"
+        assert batch_size > 0 or batch_size == -1, f"batch size should >0 or -1"
+        assert penalty in [
+            e.value for e in Penalty
+        ], f"penalty should in {[e.value for e in Penalty]}, but got {penalty}"
+        if penalty == "elasticnet":
+            assert (
+                l1_ratio > 0 and l1_ratio < 1
+            ), f"l1_ratio should >0 and <1 if use ElasticNet (L1+L2)"
+
+        self.max_iter = max_iter
+        self.eta0 = eta0
+        self.fit_intercept = fit_intercept
+        self.penalty = Penalty(penalty)
+        self.alpha = alpha
+        self.patience = patience
+        self.bsize = batch_size
+        self.early_stop = early_stop
+        if self.penalty == Penalty.L1:
+            self.l1_ratio = 1
+        elif self.penalty == Penalty.L2:
+            self.l1_ratio = 0
+        elif self.penalty == Penalty.EN:
+            self.l1_ratio = l1_ratio
+        else:
+            self.l1_ratio = None
+
+        self._w = None
+        self._b = 0.0
+
+    def _sign(self, x):
+        """
+        The sign function in perceptron.
+        if x <= 0, f(x)=-1 else if x > 0, f(x)=1
+        """
+        x = jnp.where(x <= 0, -1, 1)
+        return x
+
+    def _hinge_function(self, x, y, w, b):
+        """
+        loss = \frac{1}{n} \sum_{i=1}^{n}max\{-y_i(wx_i+b), 0\}.
+        """
+        y_hat = jnp.matmul(x, w) + b
+        y_hat = y_hat.reshape(y_hat.shape[0], 1)
+        return jnp.mean(jnp.maximum(jnp.multiply(-y_hat, y), 0))
+
+    def _update_parameters(self, x_i, y_i, w, b):
+        # if y_i * [(w * x_i) + b] <= 0, the point is misclassified used for updating w and b.
+        decision = y_i * (jnp.dot(w, x_i) + b) <= 0
+
+        dw = decision * (self.eta0 * y_i * x_i)
+        if self.fit_intercept:
+            db = decision * self.eta0 * y_i
+        else:
+            db = jnp.array([0.0])
+
+        # add regularization terms
+        if self.penalty != Penalty.NONE:
+            l1_reg = jnp.sign(w)
+            l2_reg = w
+            reg = self.l1_ratio * l1_reg + (1 - self.l1_ratio) * l2_reg
+            dw += decision * self.alpha * reg
+
+        dwb = jnp.concatenate((dw, db), axis=0)
+
+        return dwb
+
+    def fit(self, x, y):
+        """Fit Perceptron classifier.
+
+        Parameters
+        ----------
+        X : {array-like}, shape (n_samples, n_features)
+            Training data.
+
+        y : ndarray of shape (n_samples,)
+            Target values. Only -1 or 1.
+
+        Returns
+        -------
+        self : object
+            Returns an instance of self.
+        """
+
+        assert len(x.shape) == 2, f"expect x to be 2 dimension array, got {x.shape}"
+
+        n_samples, n_features = x.shape
+        if self.bsize == -1 or self.bsize > n_samples:
+            self.bsize = n_samples
+
+        w = jnp.zeros(n_features)
+        b = jnp.array([0.0])
+
+        not_early_stop = True
+
+        best_loss = jnp.inf
+        best_w = w
+        best_b = b
+        best_iter = -1
+
+        for iter in range(self.max_iter):
+            total_batch = math.ceil(float(n_samples) / self.bsize)
+            for idx in range(total_batch):
+                begin = idx * self.bsize
+                end = min((idx + 1) * self.bsize, n_samples)
+                x_slice = x[begin:end]
+                y_slice = y[begin:end]
+
+                update_w_b = jax.vmap(
+                    self._update_parameters, in_axes=(0, 0, None, None), out_axes=0
+                )(x_slice, y_slice, w, b)
+                mean_w_b = jnp.mean(jnp.array(update_w_b), axis=0)
+
+                w += not_early_stop * mean_w_b[:-1]
+
+                if self.fit_intercept:
+                    b += not_early_stop * mean_w_b[-1]
+
+            if self.early_stop:
+                sumloss = self._hinge_function(x, y, w, b)
+
+                best_loss_flag = sumloss < best_loss
+                best_loss = jnp.where(best_loss_flag, sumloss, best_loss)
+                best_w = jnp.where(best_loss_flag, w, best_w)
+                best_b = jnp.where(best_loss_flag, b, best_b)
+                best_iter = jnp.where(best_loss_flag, iter, best_iter)
+
+                not_early_stop *= iter <= best_iter + self.patience
+
+        if self.early_stop:
+            self._w = best_w
+            self._b = best_b
+        else:
+            self._w = w
+            self._b = b
+
+        return self
+
+    def predict(self, x):
+        """Perceptron estimates.
+
+        Parameters
+        ----------
+        X : {array-like}, shape (n_samples, n_features)
+            Input data for prediction.
+
+        Returns
+        -------
+        ndarray of shape (n_samples,)
+            Returns the result of the sample for each class in the model.
+        """
+        assert self._w is not None, f"the model should be trained before prediction."
+
+        pred_ = jnp.matmul(x, self._w) + self._b
+        return self._sign(pred_)


Test patch:
diff --git a/sml/linear_model/tests/BUILD.bazel b/sml/linear_model/tests/BUILD.bazel
index 2b9e7c10..7c688ab8 100644
--- a/sml/linear_model/tests/BUILD.bazel
+++ b/sml/linear_model/tests/BUILD.bazel
@@ -50,3 +50,13 @@ py_test(
         "//spu/utils:simulation",
     ],
 )
+
+py_test(
+    name = "pla_test",
+    srcs = ["pla_test.py"],
+    deps = [
+        "//sml/linear_model:pla",
+        "//spu:init",
+        "//spu/utils:simulation",
+    ],
+)
diff --git a/sml/linear_model/tests/pla_test.py b/sml/linear_model/tests/pla_test.py
new file mode 100644
index 00000000..1224fd4f
--- /dev/null
+++ b/sml/linear_model/tests/pla_test.py
@@ -0,0 +1,97 @@
+# 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
+#
+#   http://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 unittest
+
+import jax.numpy as jnp
+import pandas as pd
+import sklearn.linear_model as sk
+from sklearn.datasets import load_iris
+
+import spu.spu_pb2 as spu_pb2  # type: ignore
+import spu.utils.simulation as spsim
+from sml.linear_model.pla import Perceptron
+
+
+class UnitTests(unittest.TestCase):
+    def test_pla(self):
+        sim = spsim.Simulator.simple(
+            3, spu_pb2.ProtocolKind.ABY3, spu_pb2.FieldType.FM64
+        )
+
+        def proc(x, y):
+            model = Perceptron(
+                max_iter=20,
+                eta0=1.0,
+                penalty='elasticnet',
+                alpha=0.001,
+                fit_intercept=True,
+                l1_ratio=0.7,
+                patience=10,
+                batch_size=64,
+                early_stop=True,
+            )
+
+            return model.fit(x, y).predict(x)
+
+        def load_data():
+            iris = load_iris()
+            df = pd.DataFrame(
+                iris.data,
+                columns=['sepal length', 'sepal width', 'petal length', 'petal width'],
+            )
+            df['label'] = iris.target
+
+            # only use sepal length and sepal width features
+            # 100 samples
+            data = jnp.array(df.iloc[0:100, [0, 1, -1]])
+            x, y = data[:, :-1], data[:, -1]
+
+            # y is -1 or 1
+            y = jnp.sign(y)
+            y = jnp.where(y <= 0, -1, y)
+            y = y.reshape((y.shape[0], 1))
+
+            return x, y
+
+        # load mock data
+        x, y = load_data()
+        n_samples = len(y)
+
+        # compare with sklearn
+        sk_pla = sk.Perceptron(
+            max_iter=20,
+            eta0=1.0,
+            penalty='elasticnet',
+            alpha=0.001,
+            l1_ratio=0.7,
+            fit_intercept=True,
+        )
+        result_sk = sk_pla.fit(x, y).predict(x)
+        result_sk = result_sk.reshape(result_sk.shape[0], 1)
+        acc_sk = jnp.sum((result_sk == y)) / n_samples * 100
+
+        # run with spu
+        result = spsim.sim_jax(sim, proc)(x, y)
+        result = result.reshape(result.shape[0], 1)
+        acc_ = jnp.sum((result == y)) / n_samples * 100
+
+        # print acc
+        print(f"Accuracy in SKlearn: {acc_sk:.2f}%")
+        print(f"Accuracy in SPU: {acc_:.2f}%")
+
+
+if __name__ == "__main__":
+    unittest.main()