2013-10-19 95 views
0

我有一個矩陣D,它是m * n,我正在使用公式inv(D'* D)* D'計算僞逆,但它並不是生成的與pinv(D)相同的結果。我需要增量操作所需的術語inv(D'* D)。我的所有準確性取決於inv(D'* D),這是不正確的。有沒有其他方法可以準確地獲取inv(D'* D)?任何人都可以幫助我嗎?手動計算僞逆與matlab中的pinv不一樣

%D是我從一個博客複製的3x4矩陣,僅用於演示目的。其實原來的我也有同樣的問題,因爲它的大小太大,我不能在這裏發佈它。

D = -[1/sqrt(2) 1 1/sqrt(2) 0;0 1/sqrt(2) 1 1/sqrt(2);-1/sqrt(2) 0 1/sqrt(2) 1]; 

    B1 = pinv(D) 
    B2 = D'*inv(D*D') 

    B1 = 

     -0.353553390593274 0.000000000000000 0.353553390593274 
     -0.375000000000000 -0.176776695296637 0.125000000000000 
     -0.176776695296637 -0.250000000000000 -0.176776695296637 
     0.125000000000000 -0.176776695296637 -0.375000000000000 

    Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 
    1.904842e-017. 

    B2 = 

     -0.250000000000000     0 0.500000000000000 
     -0.500000000000000     0     0 
     0.250000000000000 -0.500000000000000     0 
         0     0 -0.750000000000000 

我需要inv(D'D)做增量操作。實際上在步驟1的問題中,每當新行被添加到D的最後一個位置時,在步驟2中,D的第一行將被刪除。所以我想使用我在這兩個步驟之前計算出的逆來找到最終的D逆。更確切地說看看這裏:

B = inv(D'*D); % if i can calculate it accurately then further work is as follows 
D1 = [D;Lr]; %Lr is last row to be added 
BLr = B-((B*Lr'*Lr*B)/(1+Lr*B*Lr')); % Row addition formula 
Fr = D1(1,:); % First row to be removed 
D2 = removerows(D1,1); 
BFr = BLr+ ((BLr*Fr'*Fr*BLr)/(1-Fr*BLr*Fr')); % row deletion formula 
B = BFr; 
Y = BFr*D2; 
+0

該警告不應忽視。考慮到惡劣條件矩陣的逆過程是一種災難。我的建議是:不惜一切代價避免使用inv()函數。這幾乎從來都不是正確的做法。 'pinv()'和反斜槓運算符通常是更好的選擇。 – nispio

+0

你使用什麼算法?我沒有時間去反編譯你的Matlab代碼。此外,第3行和第6行是否可能讀爲「'BLr = B - ((B'* Lr'* Lr * B)/(1 + Lr * B * Lr'));'」和「BFr = BLr +((BLr'* Fr'* Fr * BLr)/(1-Fr * BLr * Fr'));'「(在每個中丟失轉置?B是正方形的,所以矩陣乘法運算仍然可行,但它在當前狀態下似乎是一個不尋常的表述。 – nispio

+0

感謝#nispio。其實我試着去「B =(pinv(D'* D))* D'」,但這個也不是和pinv(D)完全一樣。另一件事是沒有轉置丟失。爲了避免反轉操作,我需要使用B進行D2反轉,這樣可以節省總體時間。在第3行中,我沒有乘以轉置,因爲在這種情況下,BLr是可以產生BFr的元素。在此之後,我做了轉置乘法Y. –

回答

4

的公式(d^TD)^ - 1 d^T或d^T(DD^T)^ - 1,你正在使用的摩爾 - 彭羅斯僞逆只如果D分別具有全列或全行排名,則爲有效。

這不是你的情況,因爲警告「矩陣接近單數」顯示。

matlab pinv命令適用於任意D,即使矩陣沒有全行或全列列。

+0

謝謝安德烈亞斯。我需要轉發一個矩陣,既不是階級赤字,但仍然存在問題。該矩陣基本上是一個矩形矩陣,我想通過(D'* D)* D'來分析某個特定目的。你能告訴我你的電子郵件ID在哪裏我可以轉發你那50 * 50大小的矩陣。 –

1

嘗試在您的矩陣上運行cond(D)並查看condition number是什麼。數字越高,矩陣越有病態。同樣,您可以運行cond(D'*D)。矩陣可以是滿秩的並且仍然是病態的。在紙面上,一個病態矩陣仍然是可逆的。但是,當您嘗試在計算機上直接反轉病態矩陣時,量化和其他效應導致的小精度錯誤可能會導致解決方案中出現難以預測的結果。

由於上述原因,通常有一種更好的方法(更穩定的數值)來實現你所要做的,而不是直接計算逆矩陣。其中許多涉及基體分解技術,如SVD。如果你能幫助我們理解爲什麼你需要inv(D'*D)這會更容易讓你指出適當的選擇方向。例如,如果您只需要僞逆,請繼續並使用pinv(),即使它與使用inv()的結果不同。 pinv()功能和​​backslash operatorinv()更具數字穩定性。

+0

條件否爲我的矩陣D出現爲5.647354603596872e + 010。爲什麼我需要inv(D'D)我已經在頂部說過了。請提出一些很好的方法來整理我的問題。 –

+0

@PrabhankarPorwal:如果D沒有像你的情況那樣滿行或列的排名,D'D或DD'是不可逆的(正如我已經說過的)。所以不管你是否需要inv(D'D)都沒有關係,它不存在。所以你的更新方法是行不通的。你必須思考爲什麼D是單數,正如nispio已經說過的那樣。 –

+0

#Andreas H.我告訴你我的D不像我的問題中指定的那樣。我知道它的(在上述問題中指定的D)等級是2,它既不是全行也不是列等級。我的實際D是50 * 50(我不能在我的問題中發帖),他的等級也是50.但條件no是5.647354603596872e + 010。 D在上面的問題中,我發佈的只是爲了展示我的問題。如果你想要我的D,我可以將它轉發給你。請提供您的電子郵件ID。 –