科列斯基分解

来自testwiki
跳转到导航 跳转到搜索

Template:NoteTA Template:Citation style 線性代數中,科列斯基分解Template:Lang-en)是指將一個正定埃爾米特矩陣分解成一個下三角矩陣與其共軛轉置乘積。這種分解方式在提高代數運算效率、蒙特卡羅方法等場合中十分有用。實數矩陣的科列斯基分解由安德烈-路易·科列斯基最先發明。實際應用中,科列斯基分解在求解線性方程組中的效率約兩倍於LU分解[1]

描述

對正定埃爾米特矩陣𝐀進行科列斯基分解,即求矩陣𝐋使下式成立

𝐀=𝐋𝐋*

其中,𝐋是一個下三角矩陣且所有對角元素均為正實數𝐋*表示𝐋的共軛轉置。每一個正定埃爾米特矩陣都有一個唯一的科列斯基分解。[2]

當矩陣𝐀是一個半正定的埃爾米特矩陣,若允許𝐋的對角線元素為零,則𝐀也存在上述形式的分解。[3]

𝐀為實數矩陣,則𝐋也為實數矩陣且科列斯基分解可改寫成𝐀=𝐋𝐋𝐓[4]

𝐀正定矩陣時,科列斯基分解是唯一的,即只存在一個對角元素均嚴格大於零的下三角矩陣,使𝐀=𝐋𝐋*成立。然而,當𝐀是半正定時,分解則不一定是唯一的。

定理的逆命題自然成立:對於某些可逆矩陣𝐋(下三角矩陣或其他矩陣),如果𝐀可被寫成𝐋𝐋*,則𝐀是一個正定的埃爾米特矩陣。

LDL分解

經典科列斯基分解的一個變形是LDL分解,即

𝐀=𝐋𝐃𝐋*

其中,𝐋是一个單位下三角矩陣𝐃是一個對角矩陣

該分解與經典科列斯基分解猶有關聯,如下:

𝐀=𝐋𝐃𝐋*=𝐋𝐃12(𝐃12)*𝐋*=𝐋𝐃12(𝐋𝐃12)*

或者,由於𝐋的對角元素必須為1,且𝐋Cholesky𝐋都是下三角矩陣,故如果已知經典科列斯基分解𝐋Cholesky,則𝐋𝐃𝐋T形式可依下式求出,設S是包含𝐋Cholesky的對角元素的對角矩陣,

𝐃=𝐒2
𝐋=𝐋Cholesky𝐒1

LDL變形如果得以有效執行,構造及使用時所需求的空間及計算的複雜性與經典科列斯基分解是相同的,但是可避免提取平方根[5]某些不存在科列斯基分解的不定矩陣,也可以執行LDL分解,此時矩陣𝐃中會出現負數元素。因此人們更傾向於使用LDL分解。對於實數矩陣,該種分解的形式可被改寫成

𝐀=𝐋𝐃𝐋𝐓

此形式通常稱為LDLT分解(或LDLT分解)。它與實對稱矩陣的特徵分解密切相關,因為對於實對稱矩陣,存在特徵分解𝐀=𝐐Λ𝐐𝐓

實例

以下乃一個實對稱矩陣的科列斯基分解︰

(41216123743164398)=(200610853)(268015003)

以下乃其LDLT分解︰

(41216123743164398)=(100310451)(400010009)(134015001)

應用

科列斯基分解主要被用於線性方程組𝐀𝐱=𝐛的求解。如果A對稱正定的,我們可以先求出𝐀=𝐋𝐋𝐓,隨後借向後替換法y求解𝐋𝐲=𝐛,再以向前替換法x求解𝐋𝐓𝐱=𝐲即得最終解。

另一種可避免在計算𝐋𝐋𝐓時需要解平方根的方法就是計算𝐀=𝐋𝐃𝐋T,然後對y求解𝐋𝐲=𝐛,最後求解𝐃𝐋T𝐱=𝐲

對於可以被改寫成對稱矩陣的線性方程組,科列斯基分解及其LDL變形是一個較高效率及較高數值穩定性的求解方法。相比之下,其效率幾近為LU分解的兩倍。Template:Fact

矩陣求逆

若欲對埃爾米特矩陣直接求逆,可以通過科列斯基分解,類似求解線性方程組一般求出逆矩陣,這需要n3次運算(12n3次乘法運算)。 該方法即便要求逐步計算也非常有效率。

非埃爾米特矩陣𝐁也可以通過下例性質求逆,因為其中𝐁𝐁*總是埃爾米特矩陣︰

𝐁1=𝐁*(𝐁𝐁*)1

計算

有各種各樣的方法用於計算科列斯基分解。 常用的演算法的計算複雜度在一般情況下為O(n3)Template:Fact

下面的算法何者較快取決於執行時的具體條件。總體而言,第一個算法會稍慢,因為其以一種不太尋常的方式讀取數據。

科列斯基算法

用於計算矩陣𝐋科列斯基算法,是以高斯消元法為基礎而調整得來的。

遞歸算法由i:=1開始,令

𝐀(1):=𝐀

在步骤i,矩陣A(i)的形式如下:

𝐀(i)=(𝐈i1000ai,i𝐛i*0𝐛i𝐁(i)),

其中,𝐈i1代表i1維的單位矩陣

此時定義矩陣𝐋i

𝐋i:=(𝐈i1000ai,i001ai,i𝐛i𝐈ni),

隨即𝐀(i)可以被改寫成

𝐀(i)=𝐋i𝐀(i+1)𝐋i*

其中

𝐀(i+1)=(𝐈i10001000𝐁(i)1ai,i𝐛i𝐛i*).

注意︰在此𝐛i𝐛i*是一個外積

重複此步驟直到i1nn步之後,我們得到A(n+1)=𝐈。因此,所求的下三角矩陣L

𝐋:=𝐋1𝐋2𝐋n.

科列斯基-巴那齊耶維茨及科列斯基-克勞特演算法

科列斯基-巴那齊耶維茨演算法以一個 5×5 矩陣執行的讀取順序(白色)及寫入順序(黄色)

寫出等式𝐀=𝐋𝐋*

𝐀=𝐋𝐋𝐓=(𝐋1100𝐋21𝐋220𝐋31𝐋32𝐋33)(𝐋11𝐋21𝐋310𝐋22𝐋3200𝐋33)=(𝐋112(symmetric)𝐋21𝐋11𝐋212+𝐋222𝐋31𝐋11𝐋31𝐋21+𝐋32𝐋22𝐋312+𝐋322+𝐋332)

則得到矩陣𝐋的元素的計算公式如下︰

𝐋j,j=Aj,jk=1j1𝐋j,k2
𝐋i,j=1𝐋j,j(Ai,jk=1j1𝐋i,k𝐋j,k),for i>j

只要𝐀是實數正定矩陣,則平方根下的表達式恆為正。

對於複埃爾米特矩陣,則適用如下公式:

𝐋j,j=Aj,jk=1j1𝐋j,k𝐋j,k*.
𝐋i,j=1𝐋j,j(𝐀i,jk=1j1𝐋i,k𝐋j,k*),for i>j.

因此,要計算𝐋i,j(ij)只需利用其的左、上方元素的值。計算通常是以以下其中一種順序進行的。

  • 科列斯基-巴那齊耶維茨演算法從矩陣L的左上角開始,依行進行計算。
  • 科列斯基-克勞特演算法從矩陣L的左上角開始,依列進行計算。

若有需要,整個矩陣可以逐個元素計算得出,無論使用何種順序讀取。

計算的穩定性

假設我們要求解一個良置的線性方程組。採用了LU分解的算法,除非進行某種繞軸旋轉,否則是不穩定的;如果算法進行了繞軸旋轉,則其誤差取決於矩陣的增長因子,這個數通常(但非總是)很小。

現在,假設算法適用科列斯基分解。如上所述,算法的效率將會是原來的兩倍。此外,無必要進行繞軸旋轉且誤差總是很小。具體而言,若要求解𝐀𝐱=𝐛𝐲表示已計算出的解,然後能解出干擾方程組(𝐀+𝐄)𝐲=𝐛,其中

𝐄2cnε𝐀2.

在這裡,2是指矩陣2-範數,而cn是一個取決於n的小常數,ε表示Template:Link-en

值得一提的是,科列斯基分解與平方根的使用有關。如果被分解的矩陣是正定的,那麽只要運算精確,矩陣中帶有平方根的元素的平方根下的數字永遠是正數。不幸的是,由於存在捨入誤差,這些數字可能為負數,並使算法擱淺。然而,此種情況僅當矩陣為病態時才會出現。一種可解決此問題的方法,是增添一個對角修正矩陣至待分解矩陣,以增加矩陣的正定性。[6] 此法雖或將減少分解的準確性,但有在某些情況卻頗有作為;譬如,執行應用於最優化的牛頓法時,若初期值距最優值較遠,則此時引入對角矩陣可以提高演算法的穩定性。

LDL分解

科列斯基分解的另一種形式——LDL分解的計算方式如下所示,

𝐀=𝐋𝐃𝐋T=(100L2110L31L321)(D1000D2000D3)(1L21L3101L32001)=(D1(symmetric)L21D1L212D1+D2L31D1L31L21D1+L32D2L312D1+L322D2+D3.)

如果 Template:Tmath 是實數矩陣,下述之遞歸計算式適用於矩陣 Template:Tmath 及 Template:Tmath 中的所有元素︰

Dj=Ajjk=1j1Ljk2Dk
Lij=1Dj(Aijk=1j1LikLjkDk),for i>j.

如果 Template:Tmath複埃爾米特矩陣,則矩陣 Template:TmathTemplate:Tmath 的計算適用以下公式:

Dj=Ajjk=1j1LjkLjk*Dk
Lij=1Dj(Aijk=1j1LikLjk*Dk),for i>j.

同樣地,若有需求,該讀取順序可以逐一計算矩陣中的每一個元素。

分塊矩陣變形

𝐀是一個不定矩陣時,LDL分解在未經過謹慎的繞軸旋轉的情況下是不穩定的;[7] 特別是,分解出的矩陣的元素可能是任意的。一个可取的改進手段是執行矩陣分塊後再執行分解,通常將原矩陣分為包含2×2的小矩陣的分塊矩陣:[8]

𝐀=𝐋𝐃𝐋T=(𝐈00𝐋21𝐈0𝐋31𝐋32𝐈)(𝐃1000𝐃2000𝐃3)(𝐈𝐋21T𝐋31T0𝐈𝐋32T00𝐈)=(𝐃1(symmetric)𝐋21𝐃1𝐋21𝐃1𝐋21T+𝐃2𝐋31𝐃1𝐋31𝐃1𝐋21T+𝐋32𝐃2𝐋31𝐃1𝐋31T+𝐋32𝐃2𝐋32T+𝐃3)

在此,矩陣中每一個元素都是一個子矩陣。故此,可依照上述遞歸公式類比計算:

𝐃j=𝐀jjk=1j1𝐋jk𝐃k𝐋jkT
𝐋ij=(𝐀ijk=1j1𝐋ik𝐃k𝐋jkT)𝐃j1

上述計算涉及矩陣相乘同精確的求逆,故而實踐中不能使用過大的分塊矩陣。

修正分解

在實踐中經常有修正科列斯基分解的需求。即,經已計算一些矩陣𝐀的科列斯基分解𝐀=𝐋𝐋*,然後在𝐀上稍作修改而產生的矩陣𝐀~,欲對其進行一個修正的科列斯基分解𝐀~=𝐋~𝐋~*。問題是,能否用已知的𝐀的分解去修正𝐀~的分解。

秩為1的修正(相加)

如果修正矩陣𝐀~=𝐀+𝐱𝐱*,則其修正的分解被稱為秩為1的修正(相加)。。

以下是一個 MATLAB 語言寫的實現秩為1的修正的小程式:

function [L] = cholupdate(L,x)
    n = length(x);
    for k=1:n
        r = sqrt(L(k,k)^2 + x(k)^2);
        c = r / L(k, k);
        s = x(k) / L(k, k);
        L(k, k) = r;
        L(k+1:n,k) = (L(k+1:n,k) + s*x(k+1:n)) / c;
        x(k+1:n) = c*x(k+1:n) - s*L(k+1:n,k);
    end
end

秩為1的修正(相減)

如果𝐀~=𝐀𝐱𝐱*,那麽只有當𝐀~仍然是正定的時候該方法才適用。 此條件下,上面的代碼也可用於相減的情況,只需要將rL(k+1:n,k)的賦值式的加號替換成减號。

證明半正定矩陣的科列斯基分解唯一

上面的演算法表明,每一個正定矩陣A都有一個科列斯基分解。此結論可以加入一些限制條件後,推廣到半正定的情況。但該條件並未被完全建立,例如,它未給出明確的數值演算法以計算科列斯基因子。

如果A是一個n×n的半正定矩陣,則序列{𝐀}={𝐀+(1/k)𝐈n}是一個由正定矩陣構成的序列。而且,在算子範數

𝐀k𝐀

在正定的情形下,每一個𝐀k都有其科列斯基分解𝐀k=𝐋k𝐋k*。根據算子範數的性質,

𝐋k2𝐋k𝐋k*=𝐀k

因而{𝐋k}是算子巴拿赫空间上的一個有界,因此是Template:Link-en(因為基礎向量空間是有限維的)。因此,它有一個帶有限制𝐋收斂子序列,也記作{𝐋k}。容易驗證,矩陣𝐋具有所需的特性,例如,滿足𝐀=𝐋𝐋*以及𝐋是下三角矩陣且其對角元素非負。對於所有的xy都有,

𝐀x,y=lim𝐀kx,y=lim𝐋k𝐋k*x,y=𝐋𝐋*x,y

因此,𝐀=𝐋𝐋*。因为基礎向量空間是有限維的,所以算子空間的所有拓撲結構都是等價的。故此,範數上𝐋k收斂於𝐋,意味着,𝐋k的每個元素都獨立地收斂於𝐋。此同時暗示,由於每個𝐋k都是對角元素非負的下三角矩陣,𝐋亦如此。

推廣

科列斯基分解可以推廣到元素中包含算子的矩陣(稱為算子矩陣)。Template:Fact{n}希爾伯特空間上的一個序列。考慮如下算子矩陣

𝐀=[𝐀11𝐀12𝐀13𝐀12*𝐀22𝐀23𝐀13*𝐀23*𝐀33]

滿足直和

=nn,

其中每一個

𝐀ij:ji

都是一個有界算子。如果𝐀是正定或半正定的,即對於有限的k及任意的

hn=1kk,

都有h,𝐀h0,則存在一個下三角的算子矩陣𝐋使得𝐀=𝐋𝐋*。此外也可以把𝐋的對角元素化為正數。

用程式語言實現

  • C語言GNU科學庫提供了幾個科列斯基分解的實現。
  • Maxima電腦代數系統︰函數Cholesky可用於計算科列斯基分解。
  • GNU Octave数值計算系統提供了一些函數以計算,修正和应用科列斯基分解。
  • LAPACK庫提供了一个高性能的科列斯基分解的實現,可以以FortranC語言及其他大多數語言讀取。
  • Python中,numpy.linalg模塊中的命令Cholesky可執行科列斯基分解。
  • Matlab中,chol命令可以簡單地對一個矩陣進行科列斯基分解。
  • R語言中,chol函數可進行科列斯基分解。
  • Wolfram Mathematica中,函數CholeskyDecomposition可以對一個矩陣執行科列斯基分解。
  • C++中,Template:Link-en中的命令chol可執行科列斯基分解。特徵庫提供了稀疏矩陣和稠密矩陣的科列斯基分解。
  • Template:Link-en中,函數Decompose提供科列斯基分解。
  • Apache Commons數學庫中有科列斯基分解的實現,可應用於JavaScala及任何其他JVM語言

參見

腳註

Template:Reflist

參考文獻

  • Template:Cite conference.
  • Template:Cite book
  • Template:Cite book.
  • S. J. Julier and J. K. Uhlmann. "A General Method for Approximating Nonlinear Transformations of ProbabilityDistributions".
  • S. J. Julier and J.K. Uhlmann, "A new extension of the Kalman filter to nonlinear systems," in Proc. AeroSense: 11th Int. Symp. Aerospace/Defence Sensing, Simulation and Controls, 1997, pp. 182–193.
  • Template:Cite book.

外部連結

科学史

  • Sur la résolution numérique des systèmes d'équations linéaires, Cholesky's 1910 manuscript, online and analyzed on BibNumTemplate:Wayback (法文) (英文) [for English, click 'A télécharger']

資訊

電腦代碼

模擬中矩陣的利用

線上計算機