pd_ipm

def pd_ipm(X, t, reg_coef, max_iter=100, tol_feas=1e-10, tol_gap=1e-6, tau_param=10, bt_params=np.array([1e-4,0.8]), display=False, bt_max_iter=1000, callback=None):
    def calc_r(lambdas, mu, v1, v2):
        return [
            t * X.dot(X.T.dot(t * lambdas)) - 1. + mu * t - v1 + v2,
            lambdas.dot(t),
            -1. / tau + v1 * lambdas,
            1. / tau + v2 * (lambdas - C)
        ]

    def update_alpha(alpha):
        return lambdas + alpha * d_lambda, mu + alpha * d_mu, v1 + alpha * d_v1, v2 + alpha * d_v2

    def backtrack(lambdas, mu, v1, v2, r_dual, r_primal, r_center1, r_center2):
        alpha = 1.
        if np.any(d_lambda > 0):
            alpha = min(alpha, min((C - lambdas[d_lambda > 0]) / d_lambda[d_lambda > 0]))
        if np.any(d_lambda < 0):
            alpha = min(alpha, min(-lambdas[d_lambda < 0] / d_lambda[d_lambda < 0]))
        if np.any(d_v1 < 0):
            alpha = min(alpha, min(-v1[d_v1 < 0] / d_v1[d_v1 < 0]))
        if np.any(d_v2 < 0):
            alpha = min(alpha, min(-v2[d_v2 < 0] / d_v2[d_v2 < 0]))
        alpha *= 0.95
        r_norm = max(norm(r_dual, np.inf), abs(r_primal), norm(r_center1, np.inf), norm(r_center2, np.inf))
        i = 0
        while i < bt_max_iter:
            i += 1
            lambdas_, mu_, v1_, v2_ = update_alpha(alpha)
            r_dual_, r_primal_, r_center1_, r_center2_ = calc_r(lambdas_, mu_, v1_, v2_)

            r_norm_ = max(norm(r_dual_, np.inf), abs(r_primal_), norm(r_center1_, np.inf), norm(r_center2_, np.inf))
            if r_norm_ < (1. - alpha * c1) * r_norm:
                return alpha
            alpha *= beta
        return alpha if r_norm_ < r_norm else None

    n, d = X.shape
    C = reg_coef / n
    c1, beta = bt_params

    lambdas = np.zeros(n)
    lambdas[t < 0] = 1. / sum(t < 0)
    lambdas[t > 0] = 1. / sum(t > 0)
    lambdas *= C / max(lambdas) / 2.
    v1 = 1. / lambdas
    v2 = 1. / (C - lambdas)
    mu = 0.
    if display:
        print ('%9s %15s %15s %15s %15s' % ('iteration', 'F', '||r_dual||', '|r_primal|', 'lambda^Tg'))

    iteration = 0
    tau = 1.
    while iteration < max_iter:
        iteration += 1
        tau = 2. * n * min(1. / tol_gap, tau_param / (v1.dot(lambdas) + v2.dot(C - lambdas)))
        r_dual, r_primal, r_center1, r_center2 = calc_r(lambdas, mu, v1, v2)
        if display:
            f = -np.sum(lambdas) + (lambdas * t).dot(X.dot(X.T.dot(lambdas * t))) / 2.
            print ('%9d %15e %15e %15e %15e' % (iteration, f, norm(r_dual, np.inf), abs(r_primal), v1.dot(lambdas) + v2.dot(C - lambdas)))
        if callback is not None:
            callback(-np.sum(lambdas) + (lambdas * t).dot(X.dot(X.T.dot(lambdas * t))) / 2.)
        if norm(r_dual, np.inf) <= tol_feas and abs(r_primal) <= tol_feas and v1.dot(lambdas) + v2.dot(C - lambdas) <= tol_gap:
            break

        Ainv = 1. / (v1 / lambdas - v2 / (lambdas - C))
        A = (X.T * t * Ainv * t).dot(X) + np.identity(d)

        y = -(r_dual + r_center1 / lambdas - r_center2 / (lambdas - C))
        y0 = r_primal

        y1 = (X.T).dot(t * Ainv * y)
        y2 = np.array(solve(A, y1)).reshape(-1)
        Pinvy = Ainv * (y - t * X.dot(y2))

        t1 = (X.T).dot(t * Ainv * t)
        t2 = np.array(solve(A, t1)).reshape(-1)
        Pinvt = Ainv * (t - t * X.dot(t2))

        d_lambda = Pinvy - Pinvt * (t.dot(Pinvy) - y0) / t.dot(Pinvt)
        d_mu = (t.dot(Pinvy) - y0) / t.dot(Pinvt)

        d_v1 = -v1 * d_lambda * 1. / lambdas + (-v1 * lambdas + 1. / tau) / lambdas
        d_v2 = -v2 * d_lambda * 1. / (lambdas - C) - (v2 * (lambdas - C) + 1. / tau) / (lambdas - C)
        alpha = backtrack(lambdas, mu, v1, v2, r_dual, r_primal, r_center1, r_center2)
        if alpha is None:
            if display:
                print("Bactracking did not found any suitable step")
            break
        lambdas, mu, v1, v2 = update_alpha(alpha)

    eps = tol_feas
    w = X.T.dot(t * lambdas)
    b = (sum([t_n - w.dot(x_n) for t_n, l_n, x_n in zip(t, lambdas, X) if l_n > eps and l_n < C - eps]) /
         sum((lambdas > eps) & (lambdas < C - eps)))
    return {'w': w, 'b': b, 'lambda': lambdas}