LASSO regression result different in Matlab and Python
Question: I am now trying to learn the ADMM algorithm (Boyd 2010) for LASSO regression.

I found out a very good example on this
[page](http://www.simonlucey.com/lasso-using-admm/).

The matlab code is shown
[here](https://dl.dropboxusercontent.com/u/22893361/code/lasso_admm.m).

I tried to convert it into python language so that I could develop a better
understanding.

Here is the code:

    
    
    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)
    

I have found the loss function did not converge after 100+ iterations. Matrix
B did not tend to be sparse, on the other hand, the matlab code worked in
different situations.

I have checked with different input data and compared with Matlab outputs, yet
I still could not get hints.

Could anybody take a try?

Thank you in advance.


Answer: My gut feeling as to why this is not working to your expectations is your
`la.solve()` call. `la.solve()` assumes that the matrix is full rank and is
independent (i.e. invertible). When you use `\` in MATLAB, what MATLAB does
under the hood is that if the matrix is full rank, the exact inverse is found.
However, should the matrix not be this way (i.e. overdetermined or
underdetermined), the solution to the system is solved by least-squares
instead. I would suggest you modify that call so that you're using
[`lstsq`](http://docs.scipy.org/doc/numpy/reference/generated/numpy.linalg.lstsq.html)
instead of `solve`. As such, simply replace your `la.solve()` call with this:

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

Note that `lstsq` returns a whole bunch of outputs in a 4-element tuple, in
addition to the solution. The solution of the system is in the first element
of this tuple, which is why I did `B = sol[0]`. What is also returned are the
sums of residues (second element), the rank (third element) and the singular
values of the matrix you are trying to invert when solving (fourth element).

* * *

Also some peculiarities that I have noticed:

  * One thing that may or may not matter is the random generation of numbers. MATLAB and Python NumPy generate random numbers differently, so this may or may not affect your solution.
  * In MATLAB, Simon Lucey's code initializes `L` to be a zero matrix such that `L = zeros(size(X));`. However, in your Python code, you initialize `L` to be this way: `L = np.zeros(C.shape);`. You are using different variables to ascertain the shape of `L`. Obviously, the code wouldn't work if there was a dimension mismatch, but that's another thing that's different. Not sure if this will affect your solution either.

* * *

So far I haven't found anything out of the ordinary, so try that fix and let
me know.

