Résultat de régression LASSO différent dans Matlab et Python

Question

J'essaie maintenant d'apprendre l'algorithme ADMM (Boyd 2010) pour la régression LASSO. J'ai découvert un très bon exemple sur ce page. Le code matlab est indiqué here. J'ai essayé de le convertir en langage Python pour que je puisse développer une meilleure compréhension.

Voici le 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) 

J'ai trouvé la fonction de perte n'a pas convergé après 100+ itérations. La matrice B n'a pas tendance à être éparse, par contre, le code matlab a fonctionné dans différentes situations.

J'ai vérifié avec différentes données d'entrée et comparé avec les sorties Matlab, mais je n'ai toujours pas pu obtenir d'indices.

Quelqu'un pourrait-il essayer?

Merci d'avance.

Answer 1

Mon intuition quant à la raison pour laquelle cela ne fonctionne pas à vos attentes est votre appel la.solve(). la.solve() suppose que la matrice est de rang entier et est indépendante (c'est-à-dire inversible). Lorsque vous utilisez \ dans MATLAB, ce que MATLAB fait sous le capot est que si la matrice est de rang complet, l'inverse exact est trouvé. Cependant, si la matrice n'est pas de cette façon (c'est-à-dire surdéterminée ou sous-déterminée), la solution au système est résolue par les moindres carrés à la place. Je vous suggère de modifier cet appel afin que vous utilisiez lstsq au lieu de solve. En tant que tel, il suffit de remplacer votre appel la.solve() avec ceci:

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

Notez que lstsq retourne un tas de sorties dans un tuple 4 éléments, en plus de la solution. La solution du système est dans le premier élément de ce tuple, c'est pourquoi j'ai fait B = sol[0]. Ce qui est également retourné, ce sont les sommes de résidus (deuxième élément), le rang (troisième élément) et les valeurs singulières de la matrice que vous essayez d'inverser lors de la résolution (quatrième élément).

---

aussi quelques particularités que je l'ai remarqué:

• Une chose qui peut ou pas d'importance est la génération aléatoire de nombres. MATLAB et Python NumPy génèrent des nombres aléatoires différemment, ce qui peut ou non affecter votre solution.
• Dans MATLAB, le code de Simon Lucey initialise L pour être une matrice zéro telle que L = zeros(size(X));. Cependant, dans votre code Python, vous initialisez L pour être ainsi: L = np.zeros(C.shape);. Vous utilisez différentes variables pour déterminer la forme de L. Évidemment, le code ne fonctionnerait pas s'il y avait une différence de dimension, mais c'est une autre chose qui est différente. Vous ne savez pas si cela affectera votre solution.

---

Jusqu'à présent, je ne l'ai pas trouvé quelque chose hors de l'ordinaire, alors essayez que fixe et laissez-moi savoir.