diff --git a/examples/portfolio/BPQP/BPQP.pdf b/examples/portfolio/BPQP/BPQP.pdf
new file mode 100644
index 0000000000..642ae9e403
Binary files /dev/null and b/examples/portfolio/BPQP/BPQP.pdf differ
diff --git a/examples/portfolio/BPQP/README.md b/examples/portfolio/BPQP/README.md
new file mode 100644
index 0000000000..12f3243205
--- /dev/null
+++ b/examples/portfolio/BPQP/README.md
@@ -0,0 +1,57 @@
+# BPQP
+The implementation of the paper: "BPQP: A Differentiable Convex Optimization Framework for Efficient End-to-End Learning" [TODO: arXiv Hyperlink]
+
+![avatar](frame.png) 
+
+# Data & Environment
+* Install python3.7, 3.8 or 3.9.
+* Install Pytorch(1.12.0 in our experiment)
+* Install the requirements in [requirements.txt](requirements.txt).
+* Install the quantitative investment platform Qlib and download the data from Qlib:
+```
+# install Qlib from source
+pip install --upgrade  cython
+git clone https://github.com/microsoft/qlib.git && cd qlib
+python setup.py install
+
+# Download the stock features of Alpha158 from Qlib
+python scripts/get_data.py qlib_data --target_dir ~/.qlib/qlib_data/cn_data --region cn --version v2 
+```
+* Run [dataset/papare_dataset.py](dataset/prepare_dataset.py) to generate train/valid/test dataset
+```
+python papare_dataset.py
+```
+# Reproduce our BPQP in large-scale QPs and LPs experiment
+
+![avatar](speed.png) 
+
+[Large scale QPs and LPs experiment.ipynb](Large scale QPs and LPs experiment.ipynb)
+
+# Reproduce our BPQP for SOCP experiment
+
+[SOCP_exp.ipynb](SOCP_exp.ipynb)
+
+# Reproduce our BPQP for end-to-end portfolio optimization
+```
+python main.py --market CN --loss e2e --predictor mlp --solver bpqp
+```
+
+# Reproduce benchmark
+* Two-Stage
+```
+python main.py --market CN --loss mse --predictor mlp --solver bpqp
+```
+
+* DC3
+```
+python main.py --market CN --loss e2e --predictor mlp --solver dc3  
+```
+
+# About the  analysis & results
+The results in the paper are located in the directory [./analysis/](./analysis/). They are named after the corresponding table or figure name in the paper.
+
+For example, if you want to obtain the results of Table 1.
+- The dependent data is located in.  [./analysis/data/Table1/](./analysis/data/Table1/)
+- The script for generating the results will be [./analysis/Table1.py](./analysis/Table1.py).
+
+Similarly, for the rest.
diff --git a/examples/portfolio/BPQP/analysis/Table1.py b/examples/portfolio/BPQP/analysis/Table1.py
new file mode 100644
index 0000000000..76b658df0c
--- /dev/null
+++ b/examples/portfolio/BPQP/analysis/Table1.py
@@ -0,0 +1,110 @@
+from pathlib import Path
+import re
+import sys
+
+import numpy as np
+import pandas as pd
+
+from data_loader import OURM, load_lp_qp_data, load_socp_data
+
+DIRNAME = Path(__file__).absolute().resolve().parent
+
+sys.path.append(str(DIRNAME))
+
+method_order = [OURM, f"{OURM}+Exact", f"{OURM}[O.F.]", "DC3", "Qpth", "QPTH", 'CVXPY', 'Exact']
+pass_order = ['F.', "B.", 'A.']
+ppname = "BPQP"
+
+
+def print_latex_table(df_mean, df_std, method_grp):
+    print("=" * 50, method_grp, "=" * 50)
+
+    def get_key(x):
+        if hasattr(x, "to_series"):  # incase of index
+            x = x.to_series()
+
+        cat_list = method_order + pass_order
+        try:
+            x = x.apply(cat_list.index)
+            return x
+        except ValueError:
+            pass
+        m = re.match(r"\d+(x|\*|×)\d+", x.iloc[0])
+        if m is not None:
+            sep = m.groups()[0]
+            x = x.apply(lambda z: int(z.split(sep)[0]))
+        return x
+
+    def rename(x):
+        # print(x)
+        if OURM in x:
+            return x.replace(OURM, ppname)
+        mp = {
+            "Qpth": "qpth/OptNet",
+            "QPTH": "qpth/OptNet",
+            "F.": "Forward",
+            "B.": "Backward",
+            "A.": "Total(Forward + Backward)",
+        }
+        return mp.get(x, x)
+
+    def tex_rename(x):
+        if 'x' in x:
+            return x.replace("x", r"×")
+        return x
+
+    def get_tex_df(df):
+        text_df2 = df
+        text_df2 = text_df2.stack(dropna=False).to_frame("time").reset_index()
+        text_df2 = text_df2.sort_values(['size', 'method', 'pass'], key=get_key)
+        tex_df = text_df2.set_index(['pass', 'method',
+                                     'size'])['time'].unstack('size').stack().unstack('pass').unstack().sort_index(
+                                         key=get_key).sort_index(axis=1, key=get_key).iloc[
+                                             ::-1,
+                                         ]
+        tex_df = tex_df.rename(rename, axis=0).rename(rename, axis=1)
+        return tex_df
+
+    def add_level(df, name, idx_name):
+        names = [idx_name] + list(df.index.names)
+        df = pd.concat({name: df})
+        df.index.names = names
+        return df
+
+    tex_df = get_tex_df(df_mean)
+    tex_df_std = get_tex_df(df_std)
+
+    scale = 10**np.ceil(np.log10(tex_df.min().min()))
+
+    def cbf(s1, s2):
+        if np.isnan(s1):
+            return "-"
+        return "{:.1f}(±{:.1f})".format(s1, s2)
+
+    def cbf2(s1, s2):
+        return s1.combine(s2, cbf)
+
+    tex_df_all = (tex_df / scale).combine((tex_df_std / scale), cbf2)
+
+    tex_df_all = tex_df_all.reindex(tex_df.index, axis=0).reindex(tex_df.columns, axis=1)
+
+    print("(scale {:.1e})".format(scale))
+
+    print(
+        add_level(add_level(tex_df_all.drop("Forward", axis=1, level=0), 'abs. time', 'metric'), method_grp,
+                  'dataset').rename(columns=tex_rename).to_latex())
+
+
+# QP or LP data
+for fname in [
+        "BPQP_QP_raw_results.csv",  # selected QP results
+        "BPQP_LP_raw_results.csv",  # Selected LP
+]:
+    method_grp = fname[5:7]
+    df_mean, df_std = load_lp_qp_data(fname)
+    df_std = df_std.reindex(df_mean.index)
+    print_latex_table(df_mean, df_std, method_grp)
+
+# ## SOCP loader
+df_mean, df_std = load_socp_data()
+print_latex_table(df_mean, df_std, "SOCP")
diff --git a/examples/portfolio/BPQP/analysis/Table2.py b/examples/portfolio/BPQP/analysis/Table2.py
new file mode 100644
index 0000000000..dfdeef59b3
--- /dev/null
+++ b/examples/portfolio/BPQP/analysis/Table2.py
@@ -0,0 +1,123 @@
+from pathlib import Path
+
+DIRNAME = Path(__file__).absolute().resolve().parent
+
+data_path = DIRNAME / "data" / "Table2"
+
+import numpy as np
+import pandas as pd
+from tqdm.auto import tqdm
+
+from e2eutils import rename, add_level, get_key, tex_rename
+
+import re
+
+
+def load_qp_acc():
+
+    acc = pd.read_csv(data_path / 'BPQP_QP_results.csv')
+
+    acc['Err.'] = acc.loc[:, "0"].apply(lambda x: re.match(r"(.*)\((.*)\)", x).groups()[0])
+    acc['Std.'] = acc.loc[:, "0"].apply(lambda x: re.match(r"(.*)\((.*)\)", x).groups()[1])
+
+    acc = acc.set_index('avg')
+
+    acc = acc.loc[acc.index.str.startswith("Acc"), :]
+
+    acc = acc.loc[~acc.index.str.startswith("Accuracy Forward"), :]
+
+    from collections import defaultdict
+
+    new_col = defaultdict(list)
+
+    for row in acc.iterrows():
+        m = re.match(r"Accuracy (?P<method>\w+) (?P<pass>\w+) ndim:(?P<var_n>\d+) neq=nineq=(?P<con_n>\d+)", row[0])
+        gd = m.groupdict()
+        gd['scale'] = f"{gd['var_n']}x{gd['con_n']}"
+        for k, v in gd.items():
+            new_col[k].append(v)
+
+    for col, values in new_col.items():
+        acc[col] = values
+
+    acc = acc.loc[:, ['scale', 'method', 'pass', 'Err.', 'Std.']]
+
+    df = acc.set_index([
+        'scale',
+        'method',
+        'pass',
+    ]).loc[:, 'Err.'].unstack('scale').swaplevel()
+    df_std = acc.set_index([
+        'scale',
+        'method',
+        'pass',
+    ]).loc[:, 'Std.'].unstack('scale').swaplevel()
+
+    # final_df = acc.rename(tex_rename, axis=1).append(dc3.swaplevel().rename(tex_rename, axis=1))
+    final_df = df.rename(tex_rename, axis=1)
+    final_df_std = df_std.rename(tex_rename, axis=1)
+
+    final_df = final_df.rename({"Backward": "QP"})
+
+    final_df_std = final_df_std.rename({"Backward": "QP"})
+    return final_df, final_df_std
+
+
+final_df, final_df_std = load_qp_acc()
+
+# load ospq
+df = pd.concat({
+    "CVXPY": pd.read_csv(data_path / 'cp_acc.csv', index_col=0),  #.iloc[50:],
+    "OSQP": pd.read_csv(data_path / 'bpqp_acc.csv', index_col=0),  # .iloc[:50],
+}, axis=1).rename(columns=rename)  # .unstack(level=0)
+
+final_df = final_df.append(pd.concat({
+    "SOCP": df.mean().unstack(),
+}, axis=0))
+
+final_df_std = final_df_std.append(pd.concat({
+    "SOCP": df.std().unstack(),
+}, axis=0))
+
+
+def show_table(final_df, final_df_std):
+
+    def agg_info(final_df):
+        final_df = final_df.astype("float")
+        final_df = final_df.sort_index(key=get_key).sort_index(key=get_key, axis=1)
+        final_df = final_df.rename(rename, axis=0)
+        final_df.columns.name = 'scale'
+        final_df = final_df.iloc[
+            ::-1,
+        ]
+
+        final_df = final_df.dropna(axis=1)
+
+        from scipy.stats import gmean
+
+        final_df = final_df.apply(gmean, axis=1).to_frame("Avg. Err.")
+        return final_df
+
+    df = agg_info(final_df)
+
+    # print(final_df.T.to_latex( float_format="{:.1e}".format, na_rep="-"))
+
+    df_std = agg_info(final_df_std)
+
+    def cbf(s1, s2):
+        if np.isnan(s1):
+            return "-"
+        return "{:.2e}(±{:.2e})".format(s1, s2)
+
+    def cbf2(s1, s2):
+        return s1.combine(s2, cbf)
+
+    df_all = df.combine(df_std, cbf2)
+
+    df_all = df_all.T.sort_index(key=get_key, axis=1)
+    df_all.columns.names = ['', "method"]
+
+    print(df_all.to_latex())
+
+
+show_table(final_df, final_df_std)
diff --git a/examples/portfolio/BPQP/analysis/Table3.py b/examples/portfolio/BPQP/analysis/Table3.py
new file mode 100644
index 0000000000..81da5ad8a3
--- /dev/null
+++ b/examples/portfolio/BPQP/analysis/Table3.py
@@ -0,0 +1,108 @@
+from pathlib import Path
+
+DIRNAME = Path(__file__).absolute().resolve().parent
+
+import numpy as np
+import pandas as pd
+from tqdm.auto import tqdm
+
+
+def read_data(path, early_stop_n=1):
+    import os
+    import pickle
+    from pathlib import Path
+    path = Path(path)
+    with open(os.path.join(path), "rb") as f:
+        df_results = pickle.load(f)
+
+    df_results = pd.DataFrame(df_results)
+
+    df_results = df_results.nlargest(early_stop_n, "valid_score")  # early stop by valid_score
+
+    focus_cols = {
+        "test_ic": "IC",
+        "test_rank_ic": "Rank IC",
+        "test_rank_icir": "Rank ICIR",
+        "test_icir": "ICIR",
+        "test_mse": "MSE",
+        "test_mdd": "MDD.",
+        'Ann.Ret.(%)': 'Ann.Ret.(%)',
+        'Sharpe': 'Sharpe',
+    }
+
+    df_results['Ann.Ret.'] = df_results['test_avg_ret'] * 240
+    df_results['Ann.Vol.'] = df_results['test_avg_std'] * 240**0.5
+    df_results['Sharpe'] = df_results['Ann.Ret.'] / df_results['Ann.Vol.']
+
+    df_results['Ann.Ret.(%)'] = df_results['Ann.Ret.'] * 100
+    return df_results[focus_cols.keys()].rename(columns=focus_cols.get)
+
+
+def read_method(key="2ST", early_stop_n=1):
+    df = []
+    for exp_path in Path(DIRNAME / f'data/Table3/{key}').glob("**/?output.dict"):
+        df.append(read_data(exp_path, early_stop_n=early_stop_n))
+    return pd.concat({key: pd.concat(df)}, axis=0)
+
+
+base_df_mlp = read_method("2ST/MLP", early_stop_n=3)
+our_df_mlp = read_method("E2E/MLP", early_stop_n=3)
+dc3_df_mlp = read_method("DC3/mlp", early_stop_n=3).append(read_method("DC3/dc3.old", early_stop_n=3)).append(
+    read_method("DC3/mlp2023-5-18/", early_stop_n=3))
+dc3_df_mlp = pd.concat({"DC3": dc3_df_mlp.droplevel(0)})
+
+selected_cols = {
+    "Prediction Metrics": ["MSE", "IC", "ICIR"],
+    "Portfolio Metrics": ["Ann.Ret.(%)", "Sharpe"],
+}
+
+
+def add_header(df):
+    return pd.concat({k: df.loc[:, cols] for k, cols in selected_cols.items()}, axis=1)
+
+
+def rename(x):
+    d = {
+        "MLP": "Two-Stage",
+        "MLP+BPQP": "COFFEE",
+        "MLP+naïve NN": "naïve NN",
+        "MLP+DC3": "DC3",
+        "Predictive Metrics": "Prediction Metrics",
+        # New:
+        "2ST": "Two-Stage",
+        "2ST/MLP": "Two-Stage",
+        "E2E/MLP": "BPQP"
+    }
+    return d.get(x, x)
+
+
+def get_key(x):
+    if isinstance(x, pd.Index):
+        x = x.to_series()
+    x = x.apply(
+        lambda i: "Two-Stage	naïve NN	DC3	BPQP	COFFEE	MSE	IC	ICIR	Ann.Ret.	Sharpe Prediction Metrics	Portfolio Metrics".
+        find(i))
+    return x
+
+
+df_all = pd.concat([add_header(base_df_mlp), add_header(our_df_mlp), add_header(dc3_df_mlp)])
+
+df_all.index.names = ['method', 'runs']
+
+from functools import partial
+
+
+def get_prec(s):
+    prec = 1 - np.floor(np.log10(s.abs().min()))
+    return prec
+
+
+def merge_runs(s, prec):
+    ft = "%%.%df" % prec
+    # print(("%s(±%s)" % (ft, ft)))
+    return ("%s(±%s)" % (ft, ft)) % (s.mean(), s.std())
+
+
+df_all = df_all.apply(lambda s: s.groupby('method').apply(partial(merge_runs, prec=get_prec(s))))
+vis_df = df_all.rename(rename).sort_index().sort_index(key=get_key, axis=1).sort_index(key=get_key)
+print(vis_df.to_latex())
diff --git a/examples/portfolio/BPQP/analysis/data/Table3/2ST/ALSTM/0output.dict b/examples/portfolio/BPQP/analysis/data/Table3/2ST/ALSTM/0output.dict
new file mode 100644
index 0000000000..a16299c465
Binary files /dev/null and b/examples/portfolio/BPQP/analysis/data/Table3/2ST/ALSTM/0output.dict differ
diff --git a/examples/portfolio/BPQP/analysis/data/Table3/2ST/ALSTM/1output.dict b/examples/portfolio/BPQP/analysis/data/Table3/2ST/ALSTM/1output.dict
new file mode 100644
index 0000000000..0d3efc4ac3
Binary files /dev/null and b/examples/portfolio/BPQP/analysis/data/Table3/2ST/ALSTM/1output.dict differ
diff --git a/examples/portfolio/BPQP/analysis/data/Table3/2ST/MLP/1/0output.dict b/examples/portfolio/BPQP/analysis/data/Table3/2ST/MLP/1/0output.dict
new file mode 100644
index 0000000000..94a1f6175a
Binary files /dev/null and b/examples/portfolio/BPQP/analysis/data/Table3/2ST/MLP/1/0output.dict differ
diff --git a/examples/portfolio/BPQP/analysis/data/Table3/2ST/MLP/1/1output.dict b/examples/portfolio/BPQP/analysis/data/Table3/2ST/MLP/1/1output.dict
new file mode 100644
index 0000000000..8135fbbe1e
Binary files /dev/null and b/examples/portfolio/BPQP/analysis/data/Table3/2ST/MLP/1/1output.dict differ
diff --git a/examples/portfolio/BPQP/analysis/data/Table3/2ST/MLP/1/2output.dict b/examples/portfolio/BPQP/analysis/data/Table3/2ST/MLP/1/2output.dict
new file mode 100644
index 0000000000..efd952b866
Binary files /dev/null and b/examples/portfolio/BPQP/analysis/data/Table3/2ST/MLP/1/2output.dict differ
diff --git a/examples/portfolio/BPQP/analysis/data/Table3/2ST/MLP/2/0output.dict b/examples/portfolio/BPQP/analysis/data/Table3/2ST/MLP/2/0output.dict
new file mode 100644
index 0000000000..abaa6207a2
Binary files /dev/null and b/examples/portfolio/BPQP/analysis/data/Table3/2ST/MLP/2/0output.dict differ
diff --git a/examples/portfolio/BPQP/analysis/data/Table3/2ST/MLP/3/0output.dict b/examples/portfolio/BPQP/analysis/data/Table3/2ST/MLP/3/0output.dict
new file mode 100644
index 0000000000..214ba08872
Binary files /dev/null and b/examples/portfolio/BPQP/analysis/data/Table3/2ST/MLP/3/0output.dict differ
diff --git a/examples/portfolio/BPQP/analysis/data/Table3/DC3/alstm/1/0output.dict b/examples/portfolio/BPQP/analysis/data/Table3/DC3/alstm/1/0output.dict
new file mode 100644
index 0000000000..b6953a1d0c
Binary files /dev/null and b/examples/portfolio/BPQP/analysis/data/Table3/DC3/alstm/1/0output.dict differ
diff --git a/examples/portfolio/BPQP/analysis/data/Table3/DC3/dc3.old/1684292312-7203455/0output.dict b/examples/portfolio/BPQP/analysis/data/Table3/DC3/dc3.old/1684292312-7203455/0output.dict
new file mode 100644
index 0000000000..5e04605599
Binary files /dev/null and b/examples/portfolio/BPQP/analysis/data/Table3/DC3/dc3.old/1684292312-7203455/0output.dict differ
diff --git a/examples/portfolio/BPQP/analysis/data/Table3/DC3/dc3.old/1684292312-7203455/1output.dict b/examples/portfolio/BPQP/analysis/data/Table3/DC3/dc3.old/1684292312-7203455/1output.dict
new file mode 100644
index 0000000000..1e0be6a99f
Binary files /dev/null and b/examples/portfolio/BPQP/analysis/data/Table3/DC3/dc3.old/1684292312-7203455/1output.dict differ
diff --git a/examples/portfolio/BPQP/analysis/data/Table3/DC3/dc3.old/1684303320-2903552/0output.dict b/examples/portfolio/BPQP/analysis/data/Table3/DC3/dc3.old/1684303320-2903552/0output.dict
new file mode 100644
index 0000000000..b3ed81bf3f
Binary files /dev/null and b/examples/portfolio/BPQP/analysis/data/Table3/DC3/dc3.old/1684303320-2903552/0output.dict differ
diff --git a/examples/portfolio/BPQP/analysis/data/Table3/DC3/dc3.old/1684303432-5532086/0output.dict b/examples/portfolio/BPQP/analysis/data/Table3/DC3/dc3.old/1684303432-5532086/0output.dict
new file mode 100644
index 0000000000..7b52957416
Binary files /dev/null and b/examples/portfolio/BPQP/analysis/data/Table3/DC3/dc3.old/1684303432-5532086/0output.dict differ
diff --git a/examples/portfolio/BPQP/analysis/data/Table3/DC3/mlp/0output.dict b/examples/portfolio/BPQP/analysis/data/Table3/DC3/mlp/0output.dict
new file mode 100644
index 0000000000..5e04605599
Binary files /dev/null and b/examples/portfolio/BPQP/analysis/data/Table3/DC3/mlp/0output.dict differ
diff --git a/examples/portfolio/BPQP/analysis/data/Table3/DC3/mlp2023-5-18/0output.dict b/examples/portfolio/BPQP/analysis/data/Table3/DC3/mlp2023-5-18/0output.dict
new file mode 100644
index 0000000000..e0dfe0208e
Binary files /dev/null and b/examples/portfolio/BPQP/analysis/data/Table3/DC3/mlp2023-5-18/0output.dict differ
diff --git a/examples/portfolio/BPQP/analysis/data/Table3/DC3/mlp2023-5-18/1output.dict b/examples/portfolio/BPQP/analysis/data/Table3/DC3/mlp2023-5-18/1output.dict
new file mode 100644
index 0000000000..334624b312
Binary files /dev/null and b/examples/portfolio/BPQP/analysis/data/Table3/DC3/mlp2023-5-18/1output.dict differ
diff --git a/examples/portfolio/BPQP/analysis/data/Table3/DC3/mlp2023-5-18/2output.dict b/examples/portfolio/BPQP/analysis/data/Table3/DC3/mlp2023-5-18/2output.dict
new file mode 100644
index 0000000000..64b64b8528
Binary files /dev/null and b/examples/portfolio/BPQP/analysis/data/Table3/DC3/mlp2023-5-18/2output.dict differ
diff --git a/examples/portfolio/BPQP/analysis/data/Table3/DC3/mlp2023-5-18/3output.dict b/examples/portfolio/BPQP/analysis/data/Table3/DC3/mlp2023-5-18/3output.dict
new file mode 100644
index 0000000000..8860fefd65
Binary files /dev/null and b/examples/portfolio/BPQP/analysis/data/Table3/DC3/mlp2023-5-18/3output.dict differ
diff --git a/examples/portfolio/BPQP/analysis/data/Table3/DC3/mlp2023-5-18/4output.dict b/examples/portfolio/BPQP/analysis/data/Table3/DC3/mlp2023-5-18/4output.dict
new file mode 100644
index 0000000000..81e8a63048
Binary files /dev/null and b/examples/portfolio/BPQP/analysis/data/Table3/DC3/mlp2023-5-18/4output.dict differ
diff --git a/examples/portfolio/BPQP/analysis/data/Table3/DC3/mlp2023-5-18/5output.dict b/examples/portfolio/BPQP/analysis/data/Table3/DC3/mlp2023-5-18/5output.dict
new file mode 100644
index 0000000000..6722ea5e65
Binary files /dev/null and b/examples/portfolio/BPQP/analysis/data/Table3/DC3/mlp2023-5-18/5output.dict differ
diff --git a/examples/portfolio/BPQP/analysis/data/Table3/E2E/ALSTM/1/0output.dict b/examples/portfolio/BPQP/analysis/data/Table3/E2E/ALSTM/1/0output.dict
new file mode 100644
index 0000000000..45171fbb34
Binary files /dev/null and b/examples/portfolio/BPQP/analysis/data/Table3/E2E/ALSTM/1/0output.dict differ
diff --git a/examples/portfolio/BPQP/analysis/data/Table3/E2E/ALSTM/1/1output.dict b/examples/portfolio/BPQP/analysis/data/Table3/E2E/ALSTM/1/1output.dict
new file mode 100644
index 0000000000..9c7e6d6290
Binary files /dev/null and b/examples/portfolio/BPQP/analysis/data/Table3/E2E/ALSTM/1/1output.dict differ
diff --git a/examples/portfolio/BPQP/analysis/data/Table3/E2E/ALSTM/1/2output.dict b/examples/portfolio/BPQP/analysis/data/Table3/E2E/ALSTM/1/2output.dict
new file mode 100644
index 0000000000..13c01c1609
Binary files /dev/null and b/examples/portfolio/BPQP/analysis/data/Table3/E2E/ALSTM/1/2output.dict differ
diff --git a/examples/portfolio/BPQP/analysis/data/Table3/E2E/ALSTM/2/0output.dict b/examples/portfolio/BPQP/analysis/data/Table3/E2E/ALSTM/2/0output.dict
new file mode 100644
index 0000000000..8872def0c9
Binary files /dev/null and b/examples/portfolio/BPQP/analysis/data/Table3/E2E/ALSTM/2/0output.dict differ
diff --git a/examples/portfolio/BPQP/analysis/data/Table3/E2E/ALSTM/2/others/0output.dict b/examples/portfolio/BPQP/analysis/data/Table3/E2E/ALSTM/2/others/0output.dict
new file mode 100644
index 0000000000..c792377717
Binary files /dev/null and b/examples/portfolio/BPQP/analysis/data/Table3/E2E/ALSTM/2/others/0output.dict differ
diff --git a/examples/portfolio/BPQP/analysis/data/Table3/E2E/MLP/1/0output.dict b/examples/portfolio/BPQP/analysis/data/Table3/E2E/MLP/1/0output.dict
new file mode 100644
index 0000000000..130ae079e1
Binary files /dev/null and b/examples/portfolio/BPQP/analysis/data/Table3/E2E/MLP/1/0output.dict differ
diff --git a/examples/portfolio/BPQP/analysis/data/Table3/E2E/MLP/1/1output.dict b/examples/portfolio/BPQP/analysis/data/Table3/E2E/MLP/1/1output.dict
new file mode 100644
index 0000000000..20a10e5e49
Binary files /dev/null and b/examples/portfolio/BPQP/analysis/data/Table3/E2E/MLP/1/1output.dict differ
diff --git a/examples/portfolio/BPQP/analysis/data/Table3/E2E/MLP/1/2output.dict b/examples/portfolio/BPQP/analysis/data/Table3/E2E/MLP/1/2output.dict
new file mode 100644
index 0000000000..af884619a9
Binary files /dev/null and b/examples/portfolio/BPQP/analysis/data/Table3/E2E/MLP/1/2output.dict differ
diff --git a/examples/portfolio/BPQP/analysis/data/Table3/E2E/MLP/2/0output.dict b/examples/portfolio/BPQP/analysis/data/Table3/E2E/MLP/2/0output.dict
new file mode 100644
index 0000000000..fba2613ba0
Binary files /dev/null and b/examples/portfolio/BPQP/analysis/data/Table3/E2E/MLP/2/0output.dict differ
diff --git a/examples/portfolio/BPQP/analysis/data/Table3/E2E/MLP/3/0output.dict b/examples/portfolio/BPQP/analysis/data/Table3/E2E/MLP/3/0output.dict
new file mode 100644
index 0000000000..01aef8a129
Binary files /dev/null and b/examples/portfolio/BPQP/analysis/data/Table3/E2E/MLP/3/0output.dict differ
diff --git a/examples/portfolio/BPQP/analysis/data_loader.py b/examples/portfolio/BPQP/analysis/data_loader.py
new file mode 100644
index 0000000000..b9c0d98e03
--- /dev/null
+++ b/examples/portfolio/BPQP/analysis/data_loader.py
@@ -0,0 +1,104 @@
+import numpy as np
+import pandas as pd
+from pathlib import Path
+import re
+
+DIRNAME = Path(__file__).absolute().resolve().parent
+OURM = "OSQP"
+
+# # Outlines: Loading Figure1 data
+def load_lp_qp_data(fname):
+
+    data_path = DIRNAME / "data/Table1"
+
+    df = pd.read_csv(data_path / fname, index_col=0)
+
+    df = df.loc[:, df.columns.str.lower().str.startswith("time")]
+
+    # Add only forward
+    for col in df.columns:
+        if OURM in col and "Forward" in col:
+            fcol = col.replace(OURM, f"{OURM}[O.F.]")
+            df[fcol] = df[col]
+            bcol = fcol.replace("Forward", "Backward")
+            df[bcol] = df[col.replace(OURM, "Exact").replace("Forward", "Backward")]
+
+    # Add All time col
+    for col in df.columns:
+        if "Forward" in col:
+            back_col = col.replace("Forward", "Backward")
+            df[col.replace("Forward", "All")] = df[col] + df[back_col]
+
+    def format_df(df):
+        df = df.to_frame("time")
+        from collections import defaultdict
+
+        new_col = defaultdict(list)
+
+        for row in df.iterrows():
+            m = re.match(
+                r"Time (?P<method>(\w|\[O\.F\.\])+) (?P<pass>\w+) ndim:(?P<var_n>\d+) neq=nineq=(?P<con_n>\d+)", row[0])
+            gd = m.groupdict()
+            gd['size'] = f"{gd['var_n']}x{gd['con_n']}"
+            for k, v in gd.items():
+                new_col[k].append(v)
+
+        for col, values in new_col.items():
+            df[col] = values
+
+        df['pass'] = df['pass'].apply(lambda x: f"{x[0]}.")
+
+        df = df.set_index(["method", "pass", "var_n", "con_n", "size"])
+
+        df = df.iloc[:, 0].unstack('pass')
+        return df
+
+    df_mean, df_std = format_df(df.mean()), format_df(df.std())
+
+    return df_mean, df_std
+
+
+# %% [markdown]
+# ## Outlines: SOCP data
+# DATA_FOLDER = "SOCP"
+DATA_FOLDER = "SOCPV02"
+
+
+## loading
+def format_data(df, method_name):
+    name_map = {
+        "10": (10, 5, "10x5"),
+        "50": (50, 10, "50x10"),
+        "100": (100, 20, "100x20"),
+        "500": (500, 100, "500x100"),
+    }
+    df = df.copy()
+    df.columns = pd.MultiIndex.from_tuples(((ps, *name_map[var_n]) for ps, var_n in df.columns))
+    df = pd.concat({method_name: df}, axis=1)
+    df.columns.names = "method	pass var_n	con_n	size".split()
+    return df
+
+
+def load_socp_data():
+    folder = DIRNAME / "data" / "Table1" / DATA_FOLDER
+    # SCS
+    df_b = pd.read_csv(folder / "cvxpy_backward_time.csv", index_col=0)  # TODO
+    df_f = pd.read_csv(folder / "cvxpy_forward_time.csv", index_col=0)  # TODO
+
+    # BPQP ECOS
+    df_b_o = pd.read_csv(folder / "bpqp_backward_time.csv", index_col=0)  # TODO
+    df_f_o = pd.read_csv(folder / "ecos_forward.csv", index_col=0)  # TODO
+
+    # Exact
+    df_b_ex = pd.read_csv(folder / "exact_back_time.csv", index_col=0)  # TODO
+
+    # BPQP Only Ours
+    df_cvxpy = format_data(pd.concat({'B.': df_b, "F.": df_f, "A.": df_f + df_b}, axis=1), "CVXPY")
+    df_bpqp = format_data(pd.concat({'B.': df_b_o, "F.": df_f_o, "A.": df_f_o + df_b_o}, axis=1), "OSQP")
+    df_ex = format_data(pd.concat({'B.': df_b_ex, "F.": df_b_ex.head(0), "A.": df_b_ex.head(0)}, axis=1), "Exact")
+    df_of = format_data(pd.concat({'B.': df_b_ex, "F.": df_f_o, "A.": df_f_o + df_b_ex}, axis=1), "OSQP[O.F.]")
+
+    df_all =  pd.concat([df_bpqp, df_cvxpy, df_ex, df_of], axis=1)
+    df_mean = df_all.mean().unstack('pass')
+    df_std = df_all.std().unstack('pass')
+    return df_mean, df_std
diff --git a/examples/portfolio/BPQP/analysis/e2eutils.py b/examples/portfolio/BPQP/analysis/e2eutils.py
new file mode 100644
index 0000000000..0611be7fa8
--- /dev/null
+++ b/examples/portfolio/BPQP/analysis/e2eutils.py
@@ -0,0 +1,58 @@
+import pandas as pd
+import re
+
+ourm = "OSQP"
+method_order = [ourm, f"{ourm}+Exact", f"{ourm}[O.F.]", "DC3", "Qpth","QPTH", "qpth/OptNet", 'CVXPY', 'Exact']
+pass_order = ['F.', "B.", 'A.']
+
+PPNAME = "BPQP"
+
+
+def rename(x):
+    # print(x)
+    if ourm in x:
+        return x.replace(ourm, PPNAME)
+    mp = {
+        "Qpth": "qpth/OptNet",
+        "QPTH": "qpth/OptNet",
+        "F.": "Forward",
+        "B.": "Backward",
+        "A.": "Total(Forward + Backward)",
+        # columns
+        "10": "10×5",
+        "50": "50×10",
+        "100": "100×20",
+        "500": "500×100",
+    }
+    return mp.get(x, x)
+
+
+
+def get_key(x):
+    if hasattr(x, "to_series"): # incase of index
+        x = x.to_series()
+
+    cat_list = method_order + pass_order
+    try:
+        x = x.apply(cat_list.index)
+        return x
+    except ValueError:
+        pass
+    m = re.match(r"\d+(x|\*|×)\d+", x.iloc[0])
+    if m is not None:
+        sep = m.groups()[0]
+        x = x.apply(lambda z: int(z.split(sep)[0]))
+    return x
+
+def tex_rename(x):
+    m = re.match(r"\d+(x|\*|×)\d+", x)
+    if m is not None:
+        sep = m.groups()[0]     
+        return x.replace(sep, r"×")
+    return x
+
+def add_level(df, name, idx_name=None):
+    names = [idx_name] + list(df.index.names)
+    df = pd.concat({name:df})
+    df.index.names = names
+    return df
diff --git a/examples/portfolio/BPQP/dataloader.py b/examples/portfolio/BPQP/dataloader.py
new file mode 100644
index 0000000000..72659a01e4
--- /dev/null
+++ b/examples/portfolio/BPQP/dataloader.py
@@ -0,0 +1,84 @@
+import torch
+import numpy as np
+import pandas as pd
+
+torch.set_default_dtype(torch.float32)
+
+
+def get_daily_code(df):
+    return df.reset_index(level=0).index.values
+
+
+class DataLoader:
+    def __init__(
+        self,
+        df_feature,
+        df_label,
+        region="CN",
+        suffix="Train",
+        batch_size=256,
+        shuffle=False,
+        device=None,
+        risk_root="./dataset/riskdata",
+    ):
+        assert len(df_feature) == len(df_label)
+        self.device = device
+
+        self.df_feature = df_feature
+        self.df_label = df_label
+
+        self.feature = torch.from_numpy(self.df_feature.values).to(self.device).to(torch.float32)
+        self.label = torch.from_numpy(self.df_label.values).to(self.device).to(torch.float32)
+
+        self.region = region
+        self.suffix = suffix
+        self.root = risk_root
+        self.index = df_label.index.to_series()
+        self.batch_size = batch_size
+        self.shuffle = shuffle
+
+        self.daily_cnt = self.index.groupby("datetime").size()
+        self.daily_index = self.index.groupby("datetime").apply(get_daily_code)
+        self.daily_date = self.daily_index.index.values
+        self.end_idx = self.daily_cnt.cumsum().astype(int)
+        self.start_idx = self.end_idx.shift(1).fillna(0).astype(int)
+
+    @property
+    def batch_length(self):
+        return len(self.label) // self.batch_size
+
+    @property
+    def daily_length(self):
+        return len(self.daily_cnt)
+
+    def iter_batch(self):
+        indices = np.arange(len(self.label))
+        if self.shuffle:
+            np.random.shuffle(indices)
+
+        for i in range(len(indices))[:: self.batch_size]:
+            if len(indices) - i < self.batch_size:
+                break
+            yield i, indices[i : i + self.batch_size]
+
+    def iter_daily(self):
+        indices = np.arange(len(self.daily_cnt))
+        for i in indices:
+            yield i, slice(self.start_idx[i], self.end_idx[i])
+
+    def get_daily(self, i, slc):
+        outs = self.feature[slc], self.label[slc][:, 0], self.get_variance(self.daily_date[i]), self.daily_index[i]
+        return outs + (self.index[slc],)
+
+    def get_batch(self, i, slc):
+        outs = self.feature[slc], self.label[slc][:, 0], torch.cov(self.feature[slc]), -1  # No date in iter batch
+        return outs + (self.index[slc],)
+
+    def get_variance(self, date):
+        date = pd.Timestamp(date)
+        path = self.root + self.region + self.suffix + "/" + date.strftime("%Y%m%d") + "/factor_exp.pkl"
+        df_variance = pd.read_pickle(path)
+        return torch.from_numpy(df_variance.values).to(self.device).to(torch.float32)
+
+    def get_daily_date(self, i):
+        return pd.Timestamp(self.daily_date[i])
diff --git a/examples/portfolio/BPQP/dataset/papare_dataset.py b/examples/portfolio/BPQP/dataset/papare_dataset.py
new file mode 100644
index 0000000000..4dc43aa6cc
--- /dev/null
+++ b/examples/portfolio/BPQP/dataset/papare_dataset.py
@@ -0,0 +1,146 @@
+import os
+import pickle
+
+import numpy as np
+import pandas as pd
+import qlib
+from qlib.config import REG_CN
+from qlib.data import D
+from qlib.data.dataset.handler import DataHandlerLP
+from qlib.model.riskmodel import StructuredCovEstimator
+from qlib.utils import init_instance_by_config
+
+provider_uri = "~/.qlib/qlib_data/cn_data"
+# provider_uri = "~/.qlib/qlib_data/us_data"  # target_dir
+# if not exists_qlib_data(provider_uri):
+#     from qlib.tests.data import GetData
+#     GetData().qlib_data(target_dir=provider_uri, region=REG_US)
+qlib.init(provider_uri=provider_uri, region=REG_CN)  # REG_US
+
+region = "CN"  # US
+market = "csi500"  # sp500
+data_handler_config = {
+    "start_time": "2008-01-01",
+    "end_time": "2020-08-01",
+    "fit_start_time": "2008-01-01",
+    "fit_end_time": "2014-12-31",
+    "instruments": market,
+    "learn_processors": [
+        {"class": "DropCol", "kwargs": {"col_list": ["VWAP0", "KUP", "KUP2", "HIGH0", "IMIN5"]}},
+        {"class": "RobustZScoreNorm", "kwargs": {"fields_group": "feature", "clip_outlier": True}},
+        {"class": "DropnaLabel"},
+        {"class": "CSZFillna", "kwargs": {"fields_group": "feature"}},
+    ],
+    "infer_processors": [
+        {"class": "DropCol", "kwargs": {"col_list": ["VWAP0", "KUP", "KUP2", "HIGH0", "IMIN5"]}},
+        {"class": "RobustZScoreNorm", "kwargs": {"fields_group": "feature", "clip_outlier": True}},
+        {"class": "CSZFillna", "kwargs": {"fields_group": "feature"}},
+    ],
+    "label": ["Ref($close, -1) / $close - 1"],
+}
+
+dataset_config = {
+    "class": "DatasetH",
+    "module_path": "qlib.data.dataset",
+    "kwargs": {
+        "handler": {
+            "class": "Alpha158",
+            "module_path": "qlib.contrib.data.handler",
+            "kwargs": data_handler_config,
+        },
+        "segments": {
+            "train": ("2008-01-01", "2014-12-31"),
+            "valid": ("2015-01-01", "2016-12-31"),
+            "test": ("2017-01-01", "2020-08-01"),
+        },
+    },
+}
+
+
+def get_daily_code(df):
+    return df.reset_index(level=0).index.values
+
+
+def robust_z_score(df):
+    return (df - df.mean()) / (0.0001 + df.std())
+
+
+def remove_ret_lim(df):
+    loc = np.argwhere((df["label"].values.squeeze() < 0.099) & (df["label"].values.squeeze() > -0.099)).squeeze()
+    return df.iloc[loc, :].fillna(method="ffill")
+
+
+def prepare_risk_data(
+    df_index, region="CN", suffix="Train", T=240, start_time="2007-01-01", riskdata_root="./riskdata"
+):
+    riskmodel = StructuredCovEstimator()
+    codes = df_index.groupby("datetime").apply(get_daily_code)
+    ret_date = codes.index.values
+    price_all = (
+        D.features(D.instruments("all"), ["$close"], start_time=start_time, end_time=ret_date[-1])
+        .squeeze()
+        .unstack(level="instrument")
+    )
+    cur_idx = np.argwhere(price_all.index == ret_date[0])[0][0]
+    assert cur_idx - T + 1 >= 0
+    for i in range(len(ret_date)):
+        date = pd.Timestamp(ret_date[i])
+        print(date)
+        ref_date = price_all.index[i + cur_idx - T + 1]
+        code = codes[i]
+        price = price_all.loc[ref_date:date, code]
+        ret = price.pct_change()
+        ret.clip(ret.quantile(0.025), ret.quantile(0.975), axis=1, inplace=True)
+
+        if suffix == "Train":
+            ret = ret.groupby("datetime").apply(robust_z_score)
+
+        try:
+            cov_estimated = riskmodel.predict(ret, is_price=False, return_decomposed_components=False)
+        except:
+            # In rare cases svd is singular, using history cov.
+            print("wrong")
+            cov_estimated = ret.cov()
+        root = riskdata_root + region + suffix + "/" + date.strftime("%Y%m%d")
+        os.makedirs(root, exist_ok=True)
+        cov_estimated.to_pickle(root + "/factor_exp.pkl")
+
+
+dataset = init_instance_by_config(dataset_config)
+df_train = dataset.prepare("train", col_set=["feature", "label"], data_key=DataHandlerLP.DK_L)
+df_valid = dataset.prepare("valid", col_set=["feature", "label"], data_key=DataHandlerLP.DK_I)
+df_test = dataset.prepare("test", col_set=["feature", "label"], data_key=DataHandlerLP.DK_I)
+
+if region == "CN":
+    df_train = remove_ret_lim(df_train)
+    df_valid = remove_ret_lim(df_valid)
+    df_test = remove_ret_lim(df_test)
+elif region == "US":
+    df_train["label"].clip(df_train["label"].quantile(0.025), df_train["label"].quantile(0.975), axis=1, inplace=True)
+    df_train = df_train.fillna(method="ffill")
+    df_valid = df_valid.fillna(method="ffill")
+    df_test = df_test.fillna(method="ffill")
+else:
+    raise NotImplementedError
+
+# train label cross-sectional z-score
+df_train["label"] = df_train["label"].groupby("datetime").apply(robust_z_score)
+
+with open(
+    "./{}_feature_dataset_market_{}_{}_start{}_end{}".format(region, market, "train", "2008-01-01", "2014-12-31"), "wb"
+) as f:
+    pickle.dump(df_train, f)
+with open(
+    "./{}_feature_dataset_market_{}_{}_start{}_end{}".format(region, market, "valid", "2015-01-01", "2016-12-31"), "wb"
+) as f:
+    pickle.dump(df_valid, f)
+with open(
+    "./{}_feature_dataset_market_{}_{}_start{}_end{}".format(region, market, "test", "2017-01-01", "2020-08-01"), "wb"
+) as f:
+    pickle.dump(df_test, f)
+print("Preparing features done!")
+
+prepare_risk_data(df_train, region=region, suffix="Train", T=240, start_time="2007-01-01")
+prepare_risk_data(df_valid, region=region, suffix="Valid", T=240, start_time="2014-01-01")
+prepare_risk_data(df_test, region=region, suffix="Test", T=240, start_time="2015-01-01")
+print("preparing risk data done!")
diff --git a/examples/portfolio/BPQP/dataset/prepare_dataset.py b/examples/portfolio/BPQP/dataset/prepare_dataset.py
new file mode 100644
index 0000000000..9349d8a418
--- /dev/null
+++ b/examples/portfolio/BPQP/dataset/prepare_dataset.py
@@ -0,0 +1,144 @@
+import os
+import pickle
+
+import numpy as np
+import pandas as pd
+import qlib
+from qlib.config import REG_CN, REG_US
+from qlib.data import D
+from qlib.data.dataset.handler import DataHandlerLP
+from qlib.model.riskmodel import StructuredCovEstimator
+from qlib.utils import init_instance_by_config
+
+provider_uri = "~/.qlib/qlib_data/cn_data"
+# provider_uri = "~/.qlib/qlib_data/us_data"  # target_dir
+# if not exists_qlib_data(provider_uri):
+#     from qlib.tests.data import GetData
+#     GetData().qlib_data(target_dir=provider_uri, region=REG_US)
+qlib.init(provider_uri=provider_uri, region=REG_CN)  # REG_US
+
+region = "CN"  # US
+market = "csi500"  # sp500
+data_handler_config = {
+    "start_time": "2008-01-01",
+    "end_time": "2020-08-01",
+    "fit_start_time": "2008-01-01",
+    "fit_end_time": "2014-12-31",
+    "instruments": market,
+    "learn_processors": [
+        {"class": "DropCol", "kwargs": {"col_list": ["VWAP0", "KUP", "KUP2", "HIGH0", "IMIN5"]}},
+        {"class": "RobustZScoreNorm", "kwargs": {"fields_group": "feature", "clip_outlier": True}},
+        {"class": "DropnaLabel"},
+        {"class": "CSZFillna", "kwargs": {"fields_group": "feature"}},
+    ],
+    "infer_processors": [
+        {"class": "DropCol", "kwargs": {"col_list": ["VWAP0", "KUP", "KUP2", "HIGH0", "IMIN5"]}},
+        {"class": "RobustZScoreNorm", "kwargs": {"fields_group": "feature", "clip_outlier": True}},
+        {"class": "CSZFillna", "kwargs": {"fields_group": "feature"}},
+    ],
+    "label": ["Ref($close, -2)/Ref($close, -1) - 1"],
+}
+
+dataset_config = {
+    "class": "DatasetH",
+    "module_path": "qlib.data.dataset",
+    "kwargs": {
+        "handler": {
+            "class": "Alpha158",
+            "module_path": "qlib.contrib.data.handler",
+            "kwargs": data_handler_config,
+        },
+        "segments": {
+            "train": ("2008-01-01", "2014-12-31"),
+            "valid": ("2015-01-01", "2016-12-31"),
+            "test": ("2017-01-01", "2020-08-01"),
+        },
+    },
+}
+
+
+def get_daily_code(df):
+    return df.reset_index(level=0).index.values
+
+
+def robust_z_score(df):
+    return (df - np.nanmean(df)) / (0.0001 + np.nanstd(df))
+
+
+def remove_ret_lim(df):
+    loc = np.argwhere((df["label"].values.squeeze() < 0.099) & (df["label"].values.squeeze() > -0.099)).squeeze()
+
+    return df.iloc[loc, :].fillna(method="ffill")
+
+
+def prepare_risk_data(
+    df_index, region="CN", suffix="Train", T=240, start_time="2007-01-01", riskdata_root="./riskdata"
+):
+    riskmodel = StructuredCovEstimator(scale_return = False)
+    codes = df_index.groupby("datetime").apply(get_daily_code)
+    ret_date = codes.index.values
+    price_all = (
+        D.features(D.instruments("all"), ["$close"], start_time=start_time, end_time=ret_date[-1])
+        .squeeze()
+        .unstack(level="instrument")
+    )
+    cur_idx = np.argwhere(price_all.index == ret_date[0])[0][0]
+    assert cur_idx - T + 1 >= 0
+    for i in range(len(ret_date)):
+        date = pd.Timestamp(ret_date[i])
+        print(date)
+        ref_date = price_all.index[i + cur_idx - T + 1]
+        code = codes[i]
+        price = price_all.loc[ref_date:date, code]
+        ret = price.pct_change().dropna(how='all')
+        ret.clip(ret.quantile(0.025), ret.quantile(0.975), axis=1, inplace=True)
+
+        ret = ret.groupby("datetime").apply(robust_z_score)
+        try:
+            cov_estimated = riskmodel.predict(ret, is_price=False, return_decomposed_components=False)
+        except ValueError:
+            print('Extreme situations: zero or one tradeable stock')
+            cov_estimated = ret.cov()
+        root = riskdata_root + region + suffix + "/" + date.strftime("%Y%m%d")
+        os.makedirs(root, exist_ok=True)
+        cov_estimated.to_pickle(root + "/factor_exp.pkl")
+
+
+dataset = init_instance_by_config(dataset_config)
+df_train = dataset.prepare("train", col_set=["feature", "label"], data_key=DataHandlerLP.DK_L)
+df_valid = dataset.prepare("valid", col_set=["feature", "label"], data_key=DataHandlerLP.DK_I)
+df_test = dataset.prepare("test", col_set=["feature", "label"], data_key=DataHandlerLP.DK_I)
+
+if region == "CN":
+    df_train = remove_ret_lim(df_train)
+    df_valid = remove_ret_lim(df_valid)
+    df_test = remove_ret_lim(df_test)
+elif region == "US":
+    df_train["label"].clip(df_train["label"].quantile(0.025), df_train["label"].quantile(0.975), axis=1, inplace=True)
+    df_train = df_train.fillna(method="ffill")
+    df_valid = df_valid.fillna(method="ffill")
+    df_test = df_test.fillna(method="ffill")
+else:
+    raise NotImplementedError
+
+# train label cross-sectional z-score
+df_train["label"] = df_train["label"].groupby("datetime").apply(robust_z_score)
+
+with open(
+    "./{}_feature_dataset_market_{}_{}_start{}_end{}".format(region, market, "train", "2008-01-01", "2014-12-31"), "wb"
+) as f:
+    pickle.dump(df_train, f)
+with open(
+    "./{}_feature_dataset_market_{}_{}_start{}_end{}".format(region, market, "valid", "2015-01-01", "2016-12-31"), "wb"
+) as f:
+    pickle.dump(df_valid, f)
+with open(
+    "./{}_feature_dataset_market_{}_{}_start{}_end{}".format(region, market, "test", "2017-01-01", "2020-08-01"), "wb"
+) as f:
+    pickle.dump(df_test, f)
+print("Preparing features done!")
+
+prepare_risk_data(df_train, region=region, suffix="Train", T=240, start_time="2007-01-01")
+prepare_risk_data(df_valid, region=region, suffix="Valid", T=240, start_time="2014-01-01")
+prepare_risk_data(df_test, region=region, suffix="Test", T=240, start_time="2015-01-01")
+print("preparing risk data done!")
diff --git a/examples/portfolio/BPQP/frame.png b/examples/portfolio/BPQP/frame.png
new file mode 100644
index 0000000000..4246a57218
Binary files /dev/null and b/examples/portfolio/BPQP/frame.png differ
diff --git a/examples/portfolio/BPQP/main.py b/examples/portfolio/BPQP/main.py
new file mode 100644
index 0000000000..855eee2d43
--- /dev/null
+++ b/examples/portfolio/BPQP/main.py
@@ -0,0 +1,343 @@
+import argparse
+import copy
+import itertools
+import os
+import pickle
+import time
+import warnings
+
+import numpy as np
+import torch
+import torch.optim as optim
+from qlib.contrib.model.pytorch_alstm import ALSTMModel
+from qlib.contrib.model.pytorch_transformer import Transformer
+from dataloader import DataLoader
+from metrics import metric_fn
+from model import MLP
+from solver import NNSolver, BPQP, QPTH, PDIPM
+from train import train_epoch, test_epoch, inference
+from utils import seed_all, dict_report, write_log
+
+torch.set_default_dtype(torch.float32)
+
+warnings.filterwarnings("ignore")
+
+device = "cuda" if torch.cuda.is_available() else "cpu"
+
+import sys
+sys.path[0] = '/home/jianming/PONet/programs/BPQP'
+
+def create_loaders(args):
+    with open(args["train_data"], "rb") as fh:
+        df_train = pickle.load(fh)
+    with open(args["valid_data"], "rb") as fh:
+        df_valid = pickle.load(fh)
+    with open(args["test_data"], "rb") as fh:
+        df_test = pickle.load(fh)
+
+    train_loader = DataLoader(
+        df_train["feature"],
+        df_train["label"],
+        region=args["market"],
+        suffix="Train",
+        shuffle=args["shuffle"],
+        device=device,
+    )
+    valid_loader = DataLoader(
+        df_valid["feature"],
+        df_valid["label"],
+        region=args["market"],
+        suffix="Valid",
+        shuffle=args["shuffle"],
+        device=device,
+    )
+    test_loader = DataLoader(
+        df_test["feature"],
+        df_test["label"],
+        region=args["market"],
+        suffix="Test",
+        shuffle=args["shuffle"],
+        device=device,
+    )
+
+    return train_loader, valid_loader, test_loader
+
+
+def create_preditor(args):
+    if args["predictor"] == "mlp":
+        predictor = MLP(d_feat=args["d_feat"], hidden_size=args["hidden_size"], num_layers=args["num_layers"])
+    elif args["predictor"] == "alstm":
+        predictor = ALSTMModel(args["d_feat"], args["hidden_size"], args["num_layers"], args["dropout"], "LSTM")
+    elif args["predictor"] == "transformer":
+        predictor = Transformer(args["d_feat"], args["hidden_size"], args["num_layers"], dropout=0.5)
+    else:
+        raise NotImplementedError
+    return predictor.to(device)
+
+
+def create_solver(args):
+    if args["solver"] == "naive_nn":
+        args["grad_step"] = False
+        solver = NNSolver(args)
+        solver = solver.to(device)
+    elif args["solver"] == "dc3":
+        args["grad_step"] = True
+        solver = NNSolver(args)
+        solver = solver.to(device)
+    elif args["solver"] == "qpth":
+        solver = QPTH
+    elif args["solver"] == "primal_dual":
+        solver = PDIPM(args)
+    else:
+        solver = BPQP(args)
+    return solver
+
+
+def train_net(train_loader, valid_loader, test_loader, metric_fn, output_path, args):
+    for times in range(args["repeat"]):
+        write_log("create preditor...")
+        model = create_preditor(args)
+        solver = create_solver(args)
+
+        if args["solver"] == "naive_nn":
+            optimizer = optim.Adam(itertools.chain(model.parameters(), solver.parameters()), lr=args["lr"])
+        else:
+            optimizer = optim.Adam(model.parameters(), lr=args["lr"])
+
+        stats = {}
+        best_score = -np.inf
+        best_epoch = 0
+        stop_round = 0
+        best_param = copy.deepcopy(model.state_dict())
+
+        for epoch in range(args["n_epochs"]):
+            write_log("Running", times, "Epoch:", epoch)
+
+            write_log("training..." + str(epoch) + "/" + str(args["n_epochs"]))
+            train_epoch(model, solver, optimizer, train_loader, args)
+            torch.save(model.state_dict(), output_path + "/model.bin.e" + str(epoch))
+            torch.save(optimizer.state_dict(), output_path + "/optimizer.bin.e" + str(epoch))
+
+            write_log("evaluating..." + str(epoch) + "/" + str(args["n_epochs"]))
+            (
+                train_loss,
+                train_regret,
+                train_mse,
+                train_score,
+                train_ic,
+                train_rank_ic,
+                train_avg_ret,
+                train_avg_std,
+                train_net_value,
+                train_mdd,
+                train_icir,
+                train_rankicir,
+            ) = test_epoch(model, solver, metric_fn, train_loader, args, prefix="Train")
+            (
+                valid_loss,
+                valid_regret,
+                valid_mse,
+                valid_score,
+                valid_ic,
+                valid_rank_ic,
+                valid_avg_ret,
+                valid_avg_std,
+                valid_net_value,
+                valid_mdd,
+                valid_icir,
+                valid_rankicir,
+            ) = test_epoch(model, solver, metric_fn, valid_loader, args, prefix="Valid")
+            (
+                test_loss,
+                test_regret,
+                test_mse,
+                test_score,
+                test_ic,
+                test_rank_ic,
+                test_avg_ret,
+                test_avg_std,
+                test_net_value,
+                test_mdd,
+                test_icir,
+                test_rankicir,
+            ) = test_epoch(model, solver, metric_fn, test_loader, args, prefix="Test")
+            write_log(
+                "epoch %d: train_loss %.6f, valid_loss %.6f, test_loss %.6f"
+                % (epoch, train_loss, valid_loss, test_loss)
+            )
+            write_log(
+                "epoch %d: train_regret %.6f, valid_regret %.6f, test_regret %.6f"
+                % (epoch, train_regret, valid_regret, test_regret)
+            )
+            write_log(
+                "epoch %d: train_mse %.6f, valid_mse %.6f, test_mse %.6f" % (epoch, train_mse, valid_mse, test_mse)
+            )
+            write_log(
+                "epoch %d: train_score %.6f, valid_score %.6f, test_score %.6f"
+                % (epoch, train_score, valid_score, test_score)
+            )
+            write_log("epoch %d: train_ic %.6f, valid_ic %.6f, test_ic %.6f" % (epoch, train_ic, valid_ic, test_ic))
+            write_log(
+                "epoch %d: train_rank_ic %.6f, valid_rank_ic %.6f, test_rank_ic %.6f"
+                % (epoch, train_rank_ic, valid_rank_ic, test_rank_ic)
+            )
+            write_log(
+                "epoch %d: train_icir %.6f, valid_icir %.6f, test_icir %.6f"
+                % (epoch, train_icir, valid_icir, test_icir)
+            )
+            write_log(
+                "epoch %d: train_rank_icir %.6f, valid_rank_icir %.6f, test_rank_icir %.6f"
+                % (epoch, train_rankicir, valid_rankicir, test_rankicir)
+            )
+            write_log(
+                "epoch %d: train_net_value %.6f, valid_net_value %.6f, test_net_value %.6f"
+                % (epoch, train_net_value, valid_net_value, test_net_value)
+            )
+            write_log(
+                "epoch %d: train_max_drawdown %.6f, valid_max_drawdown %.6f, test_max_drawdown %.6f"
+                % (epoch, train_mdd, valid_mdd, test_mdd)
+            )
+            write_log(
+                "epoch %d: train_sharpe %.6f, valid_sharpe %.6f, test_sharpe %.6f"
+                % (
+                    epoch,
+                    train_avg_ret / train_avg_std * np.sqrt(252),
+                    valid_avg_ret / valid_avg_std * np.sqrt(252),
+                    test_avg_ret / test_avg_std * np.sqrt(252),
+                )
+            )
+
+            for name in ["train", "valid", "test"]:
+                dict_report(stats, name + "_loss", eval(name + "_loss"))
+                dict_report(stats, name + "_regret", eval(name + "_regret"))
+                dict_report(stats, name + "_mse", eval(name + "_mse"))
+                dict_report(stats, name + "_score", eval(name + "_score"))
+                dict_report(stats, name + "_ic", eval(name + "_ic"))
+                dict_report(stats, name + "_rank_ic", eval(name + "_rank_ic"))
+                dict_report(stats, name + "_rank_icir", eval(name + "_rankicir"))
+                dict_report(stats, name + "_icir", eval(name + "_icir"))
+                dict_report(stats, name + "_avg_ret", eval(name + "_avg_ret"))
+                dict_report(stats, name + "_avg_std", eval(name + "_avg_std"))
+                dict_report(stats, name + "_net_value", eval(name + "_net_value"))
+                dict_report(stats, name + "_mdd", eval(name + "_mdd"))
+            # early stop
+            if valid_score > best_score:
+                best_score = valid_score
+                stop_round = 0
+                best_epoch = epoch
+                best_param = copy.deepcopy(model.state_dict())
+            else:
+                stop_round += 1
+                if stop_round >= args["early_stop"]:
+                    write_log("early stop")
+                    break
+        with open(os.path.join(output_path, str(times) + "output.dict"), "wb") as f:
+            pickle.dump(stats, f)
+
+        write_log("best score:", best_score, "@", best_epoch)
+        torch.save(best_param, output_path + "/model.bin")
+        model.load_state_dict(best_param)
+
+        infe_stats = {}
+        write_log("inference..." + str(times) + "/" + str(args["repeat"]))
+
+        pred = inference(model, solver, test_loader, args)
+        pred.to_pickle(output_path + "/pred.pkl." + "test" + str(times))
+        ic, rank_ic, avg_ret, avg_std, net_value, mdd, icir, rankicir = metric_fn(pred)
+
+        write_log("%s: IC %.6f, Rank IC %.6f, ICIR %.6f, RankICIR %.6f" % ("Test", ic, rank_ic, icir, rankicir))
+        write_log("Test", ": Sharpe ", avg_ret / avg_std * np.sqrt(252))
+        write_log("Test", ": Net_Value ", net_value)
+        write_log("Test", ": Max_Drawdown ", mdd)
+        dict_report(infe_stats, "inference_IC", ic)
+        dict_report(infe_stats, "inference_RankIC", rank_ic)
+        dict_report(infe_stats, "inference_ICIR", icir)
+        dict_report(infe_stats, "inference_RnakICIR", rankicir)
+        dict_report(infe_stats, "inference_Net_Value", net_value)
+        dict_report(infe_stats, "inference_Max_Drawdown", mdd)
+        dict_report(infe_stats, "inference_Sharpe", avg_ret / avg_std * np.sqrt(252))  # Recommend Qlib backtest
+
+        write_log("save info..." + str(times) + "/" + str(args["repeat"]))
+        with open(os.path.join(output_path, str(times) + "infe.dict"), "wb") as f:
+            pickle.dump(infe_stats, f)
+
+    write_log("finished.")
+
+
+if __name__ == "__main__":
+    parser = argparse.ArgumentParser(description="BPQP")
+    # key parameters
+    parser.add_argument("--market", type=str, default="CN", choices=["US", "CN"])
+    parser.add_argument("--loss", default="e2e", choices=["e2e", "regret", "mse", "huber_loss", "softloss"])
+    parser.add_argument("--predictor", type=str, default="mlp", choices=["mlp", "alstm", "transformer"])
+    parser.add_argument(
+        "--solver", type=str, default="bpqp", choices=["qpth", "dc3", "naive_nn", "bpqp", "primal_dual"]
+    )
+
+    # train
+    parser.add_argument("--repeat", type=int, default=1)
+    parser.add_argument("--d_feat", type=int, default=153)
+    parser.add_argument("--hidden_size", type=int, default=256)
+    parser.add_argument("--num_layers", type=int, default=2)
+    parser.add_argument("--dropout", type=float, default=0.0)
+
+    # model
+    parser.add_argument("--n_epochs", type=int, default=30)
+    parser.add_argument("--lr", type=float, default=1e-4)
+    parser.add_argument("--fre_d", type=int, default=5)
+    parser.add_argument("--gamma", type=float, default=0.1)
+    parser.add_argument("--sigma", type=float, default=1) # Risk aversion coefficient
+    parser.add_argument("--early_stop", type=int, default=5)
+    parser.add_argument("--zeta", type=int, default=4)
+
+    # data
+    parser.add_argument("--batch_size", type=int, default=256)
+    parser.add_argument("--shuffle", type=bool, default=False)
+    parser.add_argument(
+        "--train_data", default="./dataset/CN_feature_dataset_market_csi500_train_start2008-01-01_end2014-12-31"
+    )
+    parser.add_argument(
+        "--valid_data", default="./dataset/CN_feature_dataset_market_csi500_valid_start2015-01-01_end2016-12-31"
+    )
+    parser.add_argument(
+        "--test_data", default="./dataset/CN_feature_dataset_market_csi500_test_start2017-01-01_end2020-08-01"
+    )
+
+    # DC3
+    parser.add_argument("--hiddenSize", type=int, default=512)  # naive NN hidden layer
+    parser.add_argument("--max_stock", type=int, default=530)
+    parser.add_argument("--grad_step", type=bool, default=True)
+    parser.add_argument("--corrEps", type=float, default=1e-4)
+    parser.add_argument("--corrTestMaxSteps", type=int, default=10)
+    parser.add_argument("--corrMomentum", type=float, default=0)
+    parser.add_argument("--softWeightEqFrac", type=float, default=0.5)
+
+    # other
+    parser.add_argument("--seed", type=int, default=0)
+
+    args = parser.parse_args()
+    args = vars(args)
+    print(args)
+
+    seed_all(args["seed"])
+
+    save_dir = os.path.join(
+        "results",
+        str(args["market"]),
+        str(args["loss"]),
+        str(args["predictor"]),
+        str(args["solver"]),
+        str(time.time()).replace(".", "-"),
+    )
+    if not os.path.exists(save_dir):
+        os.makedirs(save_dir)
+    with open(os.path.join(save_dir, "args.dict"), "wb") as f:
+        pickle.dump(args, f)
+
+    # Load data
+    print("load data")
+    train_loader, valid_loader, test_loader = create_loaders(args)
+
+    # Run method
+    print("run training")
+    train_net(train_loader, valid_loader, test_loader, metric_fn, save_dir, args)
diff --git a/examples/portfolio/BPQP/metrics.py b/examples/portfolio/BPQP/metrics.py
new file mode 100644
index 0000000000..2c6390a2a7
--- /dev/null
+++ b/examples/portfolio/BPQP/metrics.py
@@ -0,0 +1,59 @@
+import numpy as np
+import torch
+from empyrical import max_drawdown
+
+
+def metric_fn(preds):
+    preds.index.name = "datetime"
+    preds = preds[~np.isnan(preds["label"])]
+
+    # prediction metrics
+    ic = preds.groupby("datetime").apply(lambda x: x.label.corr(x.score)).mean()
+    icir = ic / preds.groupby("datetime").apply(lambda x: x.label.corr(x.score)).std()
+    rank_ic = preds.groupby("datetime").apply(lambda x: x.label.corr(x.score, method="spearman")).mean()
+    rank_icir = rank_ic / preds.groupby("datetime").apply(lambda x: x.label.corr(x.score, method="spearman")).std()
+
+    # portfolio metrics
+    avg_ret = (
+        preds.groupby("datetime").apply(lambda x: x.label.dot(x.weight_pred) - x.label.mean()).mean()
+    )
+    avg_std = preds.groupby("datetime").apply(lambda x: x.label.dot(x.weight_pred) - x.label.mean()).std()
+    net_value = (preds.groupby("datetime").apply(lambda x: x.label.dot(x.weight_pred) - x.label.mean()) + 1).cumprod()[
+        -1
+    ]
+    ret = preds.groupby("datetime").apply(lambda x: x.label.dot(x.weight_pred) - x.label.mean())
+    mdd = max_drawdown(ret)
+
+    return ic, rank_ic, avg_ret, avg_std, net_value, mdd, icir, rank_icir
+
+
+def obj_fn(weight, rets, variance, args):
+    return weight.T @ rets - 0.5 * args["sigma"] * weight.T @ variance @ weight
+
+
+def regret_loss(weight_pred, exact_weight, pred, y, variance, args):
+    return (obj_fn(weight_pred, y, variance, args) - obj_fn(exact_weight, y, variance, args)) ** 2
+
+
+def mse_loss(pred, y, args):
+    return torch.mean((pred - y) ** 2)
+
+
+def huber_loss(weight_pred, pred, y, variance, args):
+    reg_l = mse_loss(pred, y, args) - (args['gamma']*obj_fn(weight_pred, y, variance, args))
+    if reg_l > args["zeta"] ** 2:
+        return args["zeta"] * (reg_l - args["zeta"])
+    else:
+        return reg_l
+
+
+def e2e_loss(weight_pred, exact_weight, pred, y, variance, args):
+    gamma = args["gamma"]
+    assert gamma > 0
+    return regret_loss(weight_pred, exact_weight, pred, y, variance, args)*gamma + mse_loss(pred, y, args)
+
+
+def soft_loss(weight_pred, pred, y, variance, args):
+    gamma = args["gamma"]
+    assert gamma > 0
+    return mse_loss(pred, y, args)  - (weight_pred.T @ y)*gamma
diff --git a/examples/portfolio/BPQP/model.py b/examples/portfolio/BPQP/model.py
new file mode 100644
index 0000000000..715dc02e13
--- /dev/null
+++ b/examples/portfolio/BPQP/model.py
@@ -0,0 +1,27 @@
+import torch
+from torch import nn
+
+torch.set_default_dtype(torch.float32)
+device = "cuda" if torch.cuda.is_available() else "cpu"
+
+
+class MLP(nn.Module):
+    def __init__(self, d_feat, hidden_size=128, num_layers=3, dropout=0.0):
+        super().__init__()
+
+        self.mlp = nn.Sequential()
+        for i in range(num_layers):
+            if i > 0:
+                self.mlp.add_module("drop_%d" % i, nn.Dropout(dropout))
+            self.mlp.add_module("fc_%d" % i, nn.Linear(d_feat if i == 0 else hidden_size, hidden_size))
+            self.mlp.add_module("bd_%d" % i, nn.BatchNorm1d(hidden_size))
+            self.mlp.add_module("relu_%d" % i, nn.ReLU())
+
+        self.mlp.add_module("fc_out", nn.Linear(hidden_size, 1))
+
+    def forward(self, x):
+        return self.mlp(x).squeeze()
+
+
+class ECONet:
+    pass
diff --git a/examples/portfolio/BPQP/requirements.txt b/examples/portfolio/BPQP/requirements.txt
new file mode 100644
index 0000000000..09978c907b
--- /dev/null
+++ b/examples/portfolio/BPQP/requirements.txt
@@ -0,0 +1,33 @@
+torch
+cvxpy~=1.2.1
+tqdm~=4.64.0
+pandas~=1.3.5
+audtorch~=0.6.4
+osqp~=0.6.2.post5
+scipy~=1.7.3
+qpth~=0.0.15
+PYPOWER~=5.1.15
+matplotlib~=3.5.1
+cachetools~=4.2.2
+joblib~=1.1.0
+PyYAML~=6.0
+fire~=0.4.0
+python-redis-lock~=3.7.0
+scikit-learn~=1.0.2
+requests~=2.28.1
+loguru~=0.6.0
+dill~=0.3.5.1
+redis~=4.3.4
+packaging~=21.3
+hyperopt~=0.1.2
+statsmodels~=0.13.2
+filelock~=3.7.1
+pytest~=7.1.2
+python-dateutil~=2.8.2
+lxml~=4.9.1
+beautifulsoup4~=4.11.1
+ipython~=7.34.0
+empyrical~=0.5.5
+jupyter
+seaborn
+numpy
diff --git a/examples/portfolio/BPQP/run-expers.sh b/examples/portfolio/BPQP/run-expers.sh
new file mode 100644
index 0000000000..188d7d6385
--- /dev/null
+++ b/examples/portfolio/BPQP/run-expers.sh
@@ -0,0 +1,3 @@
+nohup python main.py --market CN --loss e2e --predictor mlp --solver bpqp>> ./CNmlp_e2e_bpqp.log 2>&1 &
+nohup python main.py --market CN --loss mse --predictor mlp --solver bpqp>> ./CNmlp_mse_bpqp.log 2>&1 &
+nohup python main.py --market CN --loss e2e --predictor mlp --solver dc3>> ./CNmlp_e2e_dc3.log 2>&1 &
\ No newline at end of file
diff --git a/examples/portfolio/BPQP/solver.py b/examples/portfolio/BPQP/solver.py
new file mode 100644
index 0000000000..b70dcd7109
--- /dev/null
+++ b/examples/portfolio/BPQP/solver.py
@@ -0,0 +1,312 @@
+import operator
+import time
+from functools import reduce
+
+import numpy as np
+import osqp
+import torch
+import torch.nn.functional as F
+import qpth.qp import QPFunction
+from scipy import sparse
+from torch import nn
+from torch.autograd import Function
+
+torch.set_default_dtype(torch.double)
+device = "cuda" if torch.cuda.is_available() else "cpu"
+
+
+def osqp_interface(P, q, A, lb, ub):
+    prob = osqp.OSQP()
+    prob.setup(P, q, A, lb, ub, verbose=False, eps_abs=1e-5, eps_rel=1e-5, eps_prim_inf=1e-5, eps_dual_inf=1e-5)
+    t0 = time.time()
+    res = prob.solve()
+    return res.x, res.y, time.time() - t0
+
+
+def qp_osqp_backward(x_value, y_value, P, G, A, grad_output):
+    nineq, ndim = G.shape
+    neq = A.shape[0]
+    lambs = y_value[:nineq]  # active set
+    active_set = np.concatenate([np.argwhere(lambs > 1e-4), np.argwhere(x_value <= 1e-4)])
+    bG = G[active_set, :].squeeze()
+    bb = np.zeros(neq)
+    bh = np.zeros(len(active_set))
+    bq = -grad_output.detach().cpu().numpy()
+    osnewA = np.vstack([bG, A])
+    osnewA = sparse.csc_matrix(osnewA)
+    l_new = np.hstack([bh, bb])
+    u_new = np.hstack([bh, bb])
+    x_grad, y_grad, time_spent_backward = osqp_interface(P, bq, osnewA, l_new, u_new)
+    return x_grad, y_grad, time_spent_backward
+
+
+def create_qp_instances(P, q, G, h, A, b):
+    P, q, G, h, A, b = [x.detach().cpu().numpy() for x in [P, q, G, h, A, b]]
+    n_ineq = G.shape[0]
+    P = sparse.csc_matrix(P)
+    osA = np.vstack([G, A])
+    osA = sparse.csc_matrix(osA)
+    lb = np.hstack([np.zeros(n_ineq), 1.0]) # lower weight 0.
+    ub = np.hstack([np.ones(n_ineq), 1.0]) # upper weight 0.5
+    return P, q, osA, lb, ub
+
+
+def BPQP(args, sign=-1):
+    class BPQPmethod(Function):
+        @staticmethod
+        def forward(ctx, P, q):
+            n_dim = P.shape[0]
+            n_ineq = n_dim
+            G = torch.diag_embed(torch.ones(n_dim)).to(device)
+            h = torch.zeros(n_ineq).to(device)
+            A = torch.ones(n_dim).unsqueeze(0).to(device)
+            b = torch.tensor([1]).to(device)
+
+            _P, _q, _osA, _l, _u = create_qp_instances(P, sign * q, G, h, A, b)
+            x_value, y_value, _ = osqp_interface(_P, _q, _osA, _l, _u)
+            ctx.P = _P
+            ctx.G = G.cpu().numpy()
+            ctx.A = A.cpu().numpy()
+            yy = torch.cat(
+                [
+                    torch.from_numpy(x_value).to(device).to(torch.float32),
+                    torch.from_numpy(y_value).to(device).to(torch.float32),
+                ],
+                dim=0,
+            )
+
+            ctx.save_for_backward(yy)
+            return yy[:n_dim]
+
+        @staticmethod
+        def backward(ctx, grad_output):
+            P, G, A = ctx.P, ctx.G, ctx.A
+            ndim = P.shape[0]
+            nineq = G.shape[0]
+            yy = ctx.saved_tensors[0]
+            x_star = yy[:ndim]
+            lambda_star = yy[ndim: (ndim + nineq)]
+            x_grad, _, _ = qp_osqp_backward(
+                x_star.detach().cpu().numpy(), lambda_star.detach().cpu().numpy(), P, G, A, grad_output
+            )
+            try:
+                x_grad = torch.from_numpy(x_grad).to(torch.float32).to(device)
+            except TypeError:
+                print('No solution')
+                x_grad = None
+            grads = (None, x_grad)
+            return grads
+
+    return BPQPmethod.apply
+
+
+class CVXPY:
+    pass
+
+
+class NNSolver(nn.Module):
+    def __init__(self, args):
+        super().__init__()
+        self._args = args
+        layer_sizes = [args["max_stock"] + 1, self._args["hiddenSize"], self._args["hiddenSize"]]
+        layers = reduce(
+            operator.add,
+            [
+                [nn.Linear(a, b), nn.BatchNorm1d(b), nn.ReLU(), nn.Dropout(p=args["dropout"])]
+                for a, b in zip(layer_sizes[0:-1], layer_sizes[1:])
+            ],
+        )
+        output_dim = 1
+        layers += [nn.Linear(layer_sizes[-1], output_dim)]
+
+        self.net = nn.Sequential(*layers)
+
+    def forward(self, variance, x):
+        try:
+            x_in = torch.cat([variance, x.unsqueeze(1)], dim=1)
+        except IndexError:
+            print('Stock number < 1')
+            x_in = torch.ones(variance.shape)*(1/variance.shape[0])
+        x_in_m = torch.zeros((self._args["max_stock"], self._args["max_stock"] + 1)).to(device)
+        x_in_m[: x_in.shape[0], : x_in.shape[1]] = x_in
+        out = self.net(x_in_m).squeeze()
+        if self._args["grad_step"]:
+            G = torch.diag_embed(torch.ones(self._args["max_stock"])).to(device)
+            h = torch.zeros(self._args["max_stock"]).to(device)
+            A = torch.ones(self._args["max_stock"]).unsqueeze(0).to(device)
+            b = torch.tensor([1]).to(device)
+            out = grad_steps_all(out, A, b, G, h, self._args)
+        return F.softmax(out[: len(x)])
+
+
+def grad_steps_all(Y, A, b, G, h, args):
+    lr = args["lr"]
+    eps_converge = args["corrEps"]
+    max_steps = args["corrTestMaxSteps"]
+    momentum = args["corrMomentum"]
+    Y_new = Y
+    i = 0
+    old_Y_step = 0
+    with torch.no_grad():
+        while (
+                i == 0 or torch.max(torch.abs(A.mv(Y) - b)) > eps_converge or torch.max(G.mv(Y) - h) > eps_converge
+        ) and i < max_steps:
+            # Y_step = complete_partial(Y_new,A,b)
+            ineq_step = ineq_grad(Y_new, G, h)
+            eq_step = eq_grad(Y_new, A, b)
+            Y_step = (1 - args["softWeightEqFrac"]) * ineq_step + args["softWeightEqFrac"] * eq_step
+
+            new_Y_step = lr * Y_step + momentum * old_Y_step
+            Y_new = Y_new - new_Y_step
+
+            old_Y_step = new_Y_step
+            i += 1
+
+    return Y_new
+
+
+def complete_partial(Y, A, b):
+    rank_y = torch.linalg.matrix_rank(A)
+    Z = torch.linalg.inv(A[:, :rank_y]) @ (b - Y @ A[:, rank_y:].T).T
+    X = torch.cat([Z.T, Y], dim=1)
+    return X
+
+
+def eq_grad(Y, A, b):
+    return 2 * (A.mv(Y) - b) @ A
+
+
+def ineq_grad(Y, G, h):
+    ineq_dist = G.mv(Y) - h
+    return 2 * G.mv(ineq_dist)
+
+
+def PDIPM(args, sign=-1):
+    class QPfunction_gpuFn(Function):
+        @staticmethod
+        def forward(ctx, Q, q):
+            nz = Q.shape[0]
+            nineq = nz
+            neq = 1
+            initial_x, initial_lambs, initial_nu = (
+                (1 / nz) * torch.ones(nz).to(device),
+                torch.ones(nz).to(device),
+                torch.zeros(1).to(device),
+            )
+            _G = -torch.diag_embed(torch.ones(nz)).to(device)
+            _A = torch.ones(1, nz).to(device)
+            _h = torch.zeros(nz).to(device)
+            _b = torch.tensor(1).to(device)
+
+            x_star, lambs_star, nu_star, _, k_inverse = torch_qp_solver(
+                Q, sign * q, _G, _h, _A, _b, initial_x, initial_lambs, initial_nu, sigma=0.5, max_ite=500
+            )
+
+            ctx.nz, ctx.nineq, ctx.neq = nz, nineq, neq
+
+            ctx.save_for_backward(x_star, lambs_star, nu_star, k_inverse)
+
+            return x_star
+
+        @staticmethod
+        def backward(ctx, grad_output):
+            nz, nineq, neq = ctx.nz, ctx.nineq, ctx.neq
+            x_star, lambs_star, nu_star, k_inverse = ctx.saved_tensors
+            loss_grad = torch.cat([grad_output, torch.zeros(nineq).to(device), torch.zeros(neq).to(device)], dim=0)
+            dys = -k_inverse.mv(loss_grad)
+            dzs = dys[:nz]
+            dQs = 0.5 * (bger(dzs, x_star) + bger(x_star, dzs))
+            dqs = dzs
+            # Sparse Q,q
+            dQs[dQs < 1e-8] = 0
+            grads = (-dQs, -dqs)
+            return grads
+
+    return QPfunction_gpuFn.apply
+
+
+def solve_kkt_r(Q, q, G, h, A, b, x, lambs, nu, elips):
+    nineq = int(lambs.shape[0])
+    r_dual = q + torch.mv(Q.T, x) + torch.mv(G.T, lambs) + torch.mv(A.T, nu)
+    r_cent = torch.diag_embed(lambs).mv(torch.mv(G, x) - h).to(device) + torch.ones(nineq).to(device) * elips
+    r_prim = torch.mv(A, x) - b
+    return torch.cat([r_dual, r_cent, r_prim], dim=0)
+
+
+def solve_grad_kkt_m(Q, G, h, A, x, lambs):
+    nineq = int(G.shape[0])
+    neq = int(A.shape[0])
+    L1 = torch.cat([Q, G.T, A.T], dim=1)
+    L2 = torch.cat(
+        [torch.diag_embed(lambs).mm(G), torch.diag_embed(G.mv(x) - h), torch.zeros(nineq, neq).to(device)], dim=1
+    )
+    L3 = torch.cat([A, torch.zeros(neq, nineq).to(device), torch.zeros(neq, neq).to(device)], dim=1)
+    return torch.cat([L1, L2, L3], dim=0)
+
+
+def torch_qp_solver(Q, q, G, h, A, b, x, lambs, nu, sigma, max_ite):
+    nz = int(Q.shape[0])
+    nineq = int(G.shape[0])
+    ita = -lambs.dot((G.mv(x) - h)) / nineq
+    ita_store = []
+    ita_store.append(ita)
+    for _ in range(max_ite):
+        elips = sigma * ita
+        kkt = solve_kkt_r(Q=Q, q=q, G=G, h=h, A=A, b=b, x=x, lambs=lambs, nu=nu, elips=elips)
+        L = solve_grad_kkt_m(Q=Q, G=G, h=h, A=A, x=x, lambs=lambs)
+        # 1e-5 lower set to 0
+        k_inverse = torch.linalg.inv(L)
+        delta_y = k_inverse.mv(-kkt)  # update delta_y
+        delta_lambs = delta_y[nz: (nz + nineq)]
+        try:
+            s_max = min(1, min(-lambs[delta_lambs < 0] / delta_lambs[delta_lambs < 0]))
+        except TypeError:
+            s_max = 1
+        s = 0.99 * s_max
+        x_trail = x + s * delta_y[:nz]
+        while max(G.mv(x_trail) - h) >= 0:
+            s = 0.5 * s
+            x_trail = x + s * delta_y[:nz]
+
+        lambs_trail = lambs + s * delta_y[nz: (nz + nineq)]
+        nu_trail = nu + s * delta_y[(nz + nineq):]
+        kkt_trail = solve_kkt_r(Q=Q, q=q, G=G, h=h, A=A, b=b, x=x_trail, lambs=lambs_trail, nu=nu_trail, elips=elips)
+        while torch.norm(kkt_trail) > (1 - 0.1 * s) * torch.norm(kkt):
+            s = 0.5 * s
+            x_trail = x + s * delta_y[:nz]
+            lambs_trail = lambs + s * delta_y[nz: (nz + nineq)]
+            nu_trail = nu + s * delta_y[(nz + nineq):]
+            kkt_trail = solve_kkt_r(
+                Q=Q, q=q, G=G, h=h, A=A, b=b, x=x_trail, lambs=lambs_trail, nu=nu_trail, elips=elips
+            )  # Last KKT
+        x = x + s * delta_y[:nz]
+        lambs = lambs + s * delta_y[nz: (nz + nineq)]
+        nu = nu + s * delta_y[(nz + nineq):]
+
+        ita = -lambs.dot((G.mv(x) - h)) / nineq
+        ita_store.append(ita)
+        if (ita < 1e-3) and torch.sqrt(
+                torch.norm(delta_y[:nz], p=2) ** 2 + torch.norm(delta_y[(nz + nineq):], p=2) ** 2
+        ) < 1e-3:
+            break
+    return x, lambs, nu, ita_store, k_inverse
+
+
+def bger(x, y):
+    return x.unsqueeze(1).mm(y.unsqueeze(1).T)
+
+
+def QPTH(P, q, sign=-1):
+    # qpth tend to yield unsolved solution
+    nineq, ndim = P.shape
+    P = P + torch.eye(ndim, ndim).to(device) * 1e-6
+    G = -torch.diag_embed(torch.ones(ndim)).to(device)
+    A = torch.ones(1, ndim).to(device)
+    h = torch.zeros(ndim).to(device)
+    b = torch.tensor(1).to(device)
+    qpf = QPFunction(verbose=0, maxIter=500)
+    try:
+        qpth_x_value = qpf(P, sign * q, G, h, A, b)
+    except TypeError:
+        return F.softmax(q)
+    return qpth_x_value.squeeze()
diff --git a/examples/portfolio/BPQP/speed.png b/examples/portfolio/BPQP/speed.png
new file mode 100644
index 0000000000..030068dd2c
Binary files /dev/null and b/examples/portfolio/BPQP/speed.png differ
diff --git a/examples/portfolio/BPQP/train.py b/examples/portfolio/BPQP/train.py
new file mode 100644
index 0000000000..fbfcc97619
--- /dev/null
+++ b/examples/portfolio/BPQP/train.py
@@ -0,0 +1,162 @@
+import numpy as np
+import pandas as pd
+import torch
+from tqdm import tqdm
+from solver import BPQP
+from metrics import regret_loss, mse_loss, e2e_loss, soft_loss, huber_loss
+
+device = torch.device("cuda") if torch.cuda.is_available() else torch.device("cpu")
+
+
+def train_epoch(model, solver, optimizer, train_loader, args):
+    model.train()
+    running_loss = 0
+    for i, slc in tqdm(train_loader.iter_daily(), total=train_loader.daily_length):
+        feature, label, variance, _, _ = train_loader.get_daily(i, slc)
+
+        # predict return
+        pred = model(feature)
+
+        # differentiable solver & loss
+        if args["loss"] == "e2e":
+            weight_pred = solver(variance, pred)
+            exact_weight = solver(variance, label)
+            loss = e2e_loss(weight_pred, exact_weight, pred, label, variance, args)  # supervised by ground truth weight
+        elif args["loss"] == "regret":
+            weight_pred = solver(variance, pred)
+            exact_weight = solver(variance, label)
+            loss = regret_loss(weight_pred, exact_weight, pred, label, variance, args)
+        elif args["loss"] == "mse":
+            loss = mse_loss(pred, label, args)
+        elif args["loss"] == "huber_loss":
+            weight_pred = solver(variance, pred)
+            loss = huber_loss(weight_pred, pred, label, variance, args)
+        elif args["loss"] == "softloss":
+            weight_pred = solver(variance, pred)
+            loss = soft_loss(weight_pred, pred, label, variance, args)
+        else:
+            raise NotImplementedError
+        running_loss += loss
+
+        if i % args["fre_d"] == 0 and i > 0:
+            running_loss = running_loss / args["fre_d"]
+            running_loss.backward()
+            torch.nn.utils.clip_grad_value_(model.parameters(), 3.0)
+            optimizer.step()
+            optimizer.zero_grad()
+            running_loss = 0
+
+
+def test_epoch(model, solver, metric_fn, test_loader, args, prefix="Test"):
+    model.eval()
+
+    losses = []
+    regrets = []
+    mse = []
+    preds = []
+    if args['solver'] in ['dc3', 'naive_nn']:
+        outer_solver = BPQP(args)
+
+    for i, slc in tqdm(test_loader.iter_daily(), desc=prefix, total=test_loader.daily_length):
+        feature, label, variance, _, _ = test_loader.get_daily(i, slc)
+
+        with torch.no_grad():
+            pred = model(feature)
+
+            try:
+                if args['solver'] in ['dc3','naive_nn']:
+                    weight_pred = solver(variance, pred)
+                    exact_weight = outer_solver(variance, label)
+                else:
+                    weight_pred = solver(variance, pred)
+                    exact_weight = solver(variance, label)
+            except TypeError:
+                # Extreme situations: zero or one tradeable stock
+                print("Extreme situations: zero or one tradeable stock @", test_loader.get_daily_date(i),
+                      " Tradeable stocks:", pred.item())
+                continue
+
+            if args["loss"] == "e2e":
+                loss = e2e_loss(weight_pred, exact_weight, pred, label, variance, args)
+                regret = regret_loss(weight_pred, exact_weight, pred, label, variance, args)
+                _mse = mse_loss(pred, label, args)
+            elif args["loss"] == "regret":
+                loss = regret_loss(weight_pred, exact_weight, pred, label, variance, args)
+                regret = loss
+                _mse = mse_loss(pred, label, args)
+            elif args["loss"] == "mse":
+                # regret = regret_loss(weight_pred, exact_weight, pred, label, variance, args)
+                regret = torch.zeros(1).to(device)
+                loss = mse_loss(pred, label, args)
+                _mse = loss
+            elif args["loss"] == "huber_loss":
+                regret = regret_loss(weight_pred, exact_weight, pred, label, variance, args)
+                loss = huber_loss(weight_pred, pred, label, variance, args)
+                _mse = mse_loss(pred, label, args)
+            elif args["loss"] == "softloss":
+                regret = regret_loss(weight_pred, exact_weight, pred, label, variance, args)
+                loss = soft_loss(weight_pred, pred, label, variance, args)
+                _mse = mse_loss(pred, label, args)
+            else:
+                raise NotImplementedError
+
+            preds.append(
+                pd.DataFrame(
+                    {
+                        "score": pred.cpu().numpy(),
+                        "label": label.cpu().numpy(),
+                        "weight_pred": weight_pred.cpu().numpy(),
+                        "exact_weight": exact_weight.cpu().numpy(),
+                    },
+                    index=[test_loader.get_daily_date(i)] * len(pred),
+                )
+            )
+        regrets.append(regret.item())
+        mse.append(_mse.item())
+        losses.append(loss.item())
+    # evaluate
+    preds = pd.concat(preds, axis=0)
+    ic, rank_ic, avg_ret, avg_std, net_value, mdd, icir, rankicir = metric_fn(preds)
+
+    scores = ic
+
+
+    return (
+        np.nanmean(losses),
+        np.nanmean(regrets),
+        np.nanmean(mse),
+        scores,
+        ic,
+        rank_ic,
+        avg_ret,
+        avg_std,
+        net_value,
+        mdd,
+        icir,
+        rankicir,
+    )
+
+
+def inference(model, solver, data_loader, args):
+    model.eval()
+
+    preds = []
+    for i, slc in tqdm(data_loader.iter_daily(), total=data_loader.daily_length):
+        feature, label, variance, stock_index, _ = data_loader.get_daily(i, slc)
+        with torch.no_grad():
+            pred = model(feature)
+            weight_pred = solver(variance, pred)
+            preds.append(
+                pd.DataFrame(
+                    {
+                        "stock_index": stock_index,
+                        "score": pred.cpu().numpy(),
+                        "label": label.cpu().numpy(),
+                        "weight_pred": weight_pred.cpu().numpy(),
+                    },
+                    index=[data_loader.get_daily_date(i)] * len(pred),
+                )
+            )
+
+    preds = pd.concat(preds, axis=0)
+    return preds
diff --git a/examples/portfolio/BPQP/utils.py b/examples/portfolio/BPQP/utils.py
new file mode 100644
index 0000000000..7992eea59a
--- /dev/null
+++ b/examples/portfolio/BPQP/utils.py
@@ -0,0 +1,34 @@
+import datetime
+import logging
+import random
+
+import numpy as np
+import torch
+
+
+def seed_all(seed):
+    torch.manual_seed(seed)
+    torch.cuda.manual_seed_all(seed)
+    np.random.seed(seed)
+    random.seed(seed)
+    torch.backends.cudnn.deterministic = True
+
+
+def dict_report(stats, key, value, op="append"):
+    if key in stats.keys():
+        if op == "append":
+            stats[key] = np.append(stats[key], value)
+        elif op == "concat":
+            stats[key] = np.concatenate((stats[key], value), axis=0)
+        else:
+            raise NotImplementedError
+    else:
+        stats[key] = value
+
+
+def write_log(*info):
+    # print with UTC+8 time
+    time = "[" + str(datetime.datetime.utcnow() + datetime.timedelta(hours=8))[:19] + "] -"
+    print(time, *info)
+
+    logging.info("{time} {info}".format(time=time,info=info))