0

我現在正在嘗試學習LASSO迴歸的ADMM算法(Boyd 2010)。Matlab和Python中LASSO迴歸結果不同

我在page上發現了一個很好的例子。

該matlab代碼顯示爲here

我試圖將它轉換成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輸出進行比較,但我仍然無法獲得提示。

有人可以試試嗎?

預先感謝您。

+0

請通過添加*完整*示例來更好地解決這個問題,特別是包括傳遞給函數的X,A和伽馬值。你說你使用了各種各樣的,但至少提供一套,以便其他人可以快速檢查你的代碼。 (-1不是我) – Unapiedra 2014-09-03 16:13:09

+0

感謝您的評論。我用兩個輸入文件來測試上面的代碼。請參閱[C.mat](https://www.dropbox.com/s/g0vb3s3cib614pm/C.mat?dl=0)和[lasso.mat](https://www.dropbox.com/s/57ia207tjzp4ic6/ lasso.mat?DL = 0)。請注意,這個版本與原始Matlab代碼有些不同,後者使用隨機矩陣。 – SpencerLo 2014-09-03 18:23:25

回答

1

我的直覺爲什麼這不符合您的期望是您撥打la.solve()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返回一大堆的輸出在4元元組,除了解決方案。系統的解決方案在這個元組的第一個元素,這就是爲什麼我做了B = sol[0]。同樣返回的是殘差(第二個元素),秩(第三個元素)和您在求解時嘗試反演的矩陣的奇異值的總和(第四個元素)。


也有一些特點,我已經注意到:

  • 一件事,可能或可能不會不管是隨機生成的數字。 MATLAB和Python NumPy以不同方式生成隨機數,所以這可能會影響您的解決方案。
  • 在MATLAB中,Simon Lucey的代碼將L初始化爲零矩陣,即L = zeros(size(X));。但是,在你的Python代碼中,你初始化L就是這樣的:L = np.zeros(C.shape);。您正在使用不同的變量來確定L的形狀。顯然,如果尺寸不匹配,代碼將不起作用,但這是另一回事。不知道這是否會影響您的解決方案。

到目前爲止,我還沒有發現任何不尋常的,所以嘗試修復,讓我知道。

+0

謝謝你的熱心幫助。 1.我曾嘗試使用lstsq,但它產生的結果與之前的解決方案相同。 2.我的朋友幫我修改了Python代碼,因爲它現在需要某些輸入數據(請參閱上面的鏈接評論)。我可以看到損失函數在114次迭代之後開始發散。我的觀點是,在這種情況下,隨機分配似乎不成問題。我在上面的代碼中修改了L尺寸。 – SpencerLo 2014-09-03 18:23:57

+0

@SpencerLo - 發生了什麼?我建議做什麼工作? – rayryeng 2014-09-03 18:24:42