Matlab和Python中的LASSO回归结果不同

Question

我现在正在学习ADMM算法（Boyd 2010），用于LASSO回归。

我在这个页面上找到了一个非常好的例子。

这里有个matlab代码链接。

我尝试把它转换成python语言，这样我可以更好地理解。

以下是代码：

import scipy.io as io
import scipy.sparse as sp
import scipy.linalg as la
import numpy as np

def l1_norm(x):
    return np.sum(np.abs(x))

def l2_norm(x):
    return np.dot(x.ravel().T, x.ravel())

def fast_threshold(x, threshold):
    return np.multiply(np.sign(x), np.fmax(abs(x) - threshold, 0))

def lasso_admm(X, A, gamma):
    c = X.shape[1]
    r = A.shape[1]

    C = io.loadmat("C.mat")["C"]

    L = np.zeros(X.shape)

    rho = 1e-4
    maxIter = 200
    I = sp.eye(r)
    maxRho = 5

    cost = []

    for n in range(maxIter):
        B = la.solve(np.dot(A.T, A) + rho * I, np.dot(A.T, X) + rho * C - L)

        C = fast_threshold(B + L / rho, gamma / rho)

        L = L + rho * (B - C);  

        rho = min(maxRho, rho * 1.1); 

        cost.append(0.5 * l2_norm(X - np.dot(A, B)) + gamma * l1_norm(B))

    cost = np.array(cost).ravel()

    return B, cost

data = io.loadmat("lasso.mat")
A = data["A"]
X = data["X"]    
B, cost = lasso_admm(X, A, gamma)

我发现损失函数在100多次迭代后没有收敛。矩阵B没有变得稀疏，另一方面，matlab代码在不同情况下是有效的。

我尝试了不同的输入数据，并与Matlab的输出进行了比较，但仍然没有找到线索。

有没有人可以试一下？

提前谢谢你们。

Answer 1

我觉得你遇到的问题可能和你使用的 la.solve() 有关。这个函数假设你的矩阵是满秩的，也就是说它是独立的（也就是可以逆的）。在 MATLAB 中，当你使用 \ 时，如果矩阵是满秩的，MATLAB 会找到它的确切逆矩阵。但是，如果矩阵不是满秩的（比如说过多方程或方程太少），那么系统的解会通过最小二乘法来求解。因此，我建议你把这个调用改成使用 lstsq，而不是 solve。所以，你只需要把 la.solve() 替换成这个：

sol = la.lstsq(np.dot(A.T, A) + rho * I, np.dot(A.T, X) + rho * C - L)
B = sol[0]

需要注意的是，lstsq 会返回一堆输出，组成一个四个元素的元组，除了你想要的解。系统的解在这个元组的第一个元素里，所以我用了 B = sol[0]。另外返回的还有残差的和（第二个元素）、秩（第三个元素）以及你在求解时试图逆的矩阵的奇异值（第四个元素）。

---

我还注意到一些特别的地方：

• 一个可能会影响结果的因素是随机数的生成。MATLAB 和 Python 的 NumPy 生成随机数的方式不同，所以这可能会影响你的解，也可能不会。
• 在 MATLAB 中，Simon Lucey 的代码是这样初始化 L 的：L = zeros(size(X));。而在你的 Python 代码中，你是这样初始化 L 的：L = np.zeros(C.shape);。你使用了不同的变量来确定 L 的形状。如果维度不匹配，代码肯定会出错，但这也是一个不同之处。不确定这是否会影响你的解。

---

到目前为止，我没有发现什么异常，所以试试这个修复方法，然后告诉我结果。