Cholesky 正定矩陣分解

最近在搞 UKF，需要求解 P 矩陣的開跟號，所以找了一篇文章線代啟示錄 - Cholesky 分解參考，賺寫一個 Cholesky 矩陣分解的程式，暫時先用 MATLAB 實現，方便之後再移植到單晶片上做處理。

簡單講一下原理，假定有一個正定矩陣 A，我們希望他可以被分解成一個 P 矩陣乘上自己的轉置，而這個 P 矩陣是一個下三角矩陣，也就是我們想要求得的分解結果，如下圖所示，

$\large \begin{center} A_{n\times n}=PP^{T}\\ \end{center} \begin{center} A=\begin{bmatrix} a_{11}&a_{12}&a_{13}&\cdots&a_{1n}\\ a_{21}&a_{22}&a_{23}&\cdots&a_{2n}\\ a_{31}&a_{32}&a_{33}&\cdots&a_{3n}\\ \vdots&\vdots&\vdots&\ddots&\vdots\\ a_{n1}&a_{n2}&a_{n2}&\cdots&a_{nn}\\ \end{bmatrix}, P=\begin{bmatrix} p_{11}&0&0&\cdots&0\\ p_{21}&p_{22}&0&\cdots&0\\ p_{31}&p_{32}&p_{33}&\cdots&0\\ \vdots&\vdots&\vdots&\ddots&\vdots\\ p_{n1}&p_{n2}&p_{n3}&\cdots&p_{nn}\\ \end{bmatrix} \end{center}$

如果 A 矩陣是一個正定矩陣，則可以表示成下面形式，L 和 U 分別為下三角和上三角矩陣，X 為對角矩陣，並且 A 矩陣的轉置會等於 A 矩陣，
　
$\large \begin{align*} &A=LXU\\ &A^{T}=L^{T}X^{T}U^{T}=L^{T}XU^{T} \end{align*}$

藉由 A 矩陣的轉置會等於 A 矩陣的條件，可以得到下列關係，並且帶回前式整理，

$\large \begin{align*} &A=A^{T}\\ &\Rightarrow L=U^{T},\;U=L^{T}\\ &\Rightarrow A=LXL^{T}=(LX^{\frac{1}{2}})(X^{\frac{1}{2}}L^{T})=(LX^{\frac{1}{2}})(LX^{\frac{1}{2}})^{T}=PP^{T} \end{align*}$

可以得到 A 矩陣和 P 矩陣的關係，這裡以 3x3 的矩陣做展開，再過歸納，

$\large \begin{align*} &A_{3\times 3}=PP^{T}\\ &\Rightarrow \begin{bmatrix} a_{11}&a_{12}&a_{13}\\ a_{21}&a_{22}&a_{23}\\ a_{31}&a_{32}&a_{33}\\ \end{bmatrix}= \begin{bmatrix} p_{11}&0&0\\ p_{21}&p_{22}&0\\ p_{31}&p_{32}&p_{33}\\ \end{bmatrix}\times \begin{bmatrix} p_{11}&p_{21}&p_{31}\\ 0&p_{22}&p_{32}\\ 0&0&p_{33}\\ \end{bmatrix}\\\\ &a_{11}=p_{11}^{2}\\ &a_{21}=p_{11}\,p_{21}\\ &a_{31}=p_{11}\,p_{31}\\ &a_{22}=p_{21}^{2}+p_{22}^{2}\\ &a_{32}=p_{21}\,p_{31}+p_{22}\,p_{32}\\ &a_{33}=p_{31}^{2}+p_{32}^{2}+p_{33}^{2}\\ \end{align*}$

上式反推得到下列關係，

$\large \begin{align*} &p_{11}=(a_{11})^{1/2}\\ &p_{21}=a_{21}/p_{11}\\ &p_{31}=a_{31}/p_{11}\\ &p_{22}=(a_{22}-p_{21}^{2})^{1/2}\\ &p_{32}=(a_{32}-p_{21}\,p_{31})/p_{22}\\ &p_{33}=(a_{33}-p_{31}^{2}-p_{32}^{2})^{1/2}\\ \end{align*}$

透過上述關係可歸納出在 n x n 的正定矩陣分解出來的 P 矩陣如下。

$\large \begin{align*} &p_{ij}= \begin{cases} (a_{ij}-\sum_{k=1}^{i-1}p_{ik}^{2})^{1/2},\quad\qquad i=j\\ (a_{ij}-\sum_{k=1}^{j-1}p_{ik}\,p_{jk})/p_{jj},\quad i\neq j \end{cases}\\ &i,\,j=1,\,2,\,\cdots,\,n \end{align*}$

最後透過 MATLAB 來做驗證，MATLAB 有一個 Chol(．) 的指令，用來做 cholesky 分解的，下面自己撰寫好 cholesky 分解後，透過 A = P P' 來驗證結果，當然也可以透過 Chol(．) 來判斷是否正確，程式沒有很完善，在正定矩陣的判斷等等沒有處理，如果要實際應用，建議需要另外加入一些條件處理特殊狀況。

cholesky.m link

function sqrt_matrix = cholesky( pdMatrix )

[n, m] = size(pdMatrix);
sqrt_matrix = zeros(n);
for j = 1 : n
    for i = j : n
        if i == j
            tempSum = 0;
            for k = 1 : i - 1
                tempSum = tempSum + sqrt_matrix(i, k)^2;
            end
            sqrt_matrix(i, j) = sqrt(pdMatrix(i, j) - tempSum);
            else
                tempSum = 0;
        for k = 1 : j - 1
                    tempSum = tempSum + sqrt_matrix(i, k) * sqrt_matrix(j, k);
                end
                sqrt_matrix(i, j) = (pdMatrix(i, j) - tempSum) / sqrt_matrix(j, j);
            end
        end
    end

end

testCholesky.m link

n = 10;
X = diag(rand(n, 1));
U = orth(rand(n, n));

pdMatrix = U' * X * U;

P_U = chol(pdMatrix, 'upper');
P_L = chol(pdMatrix, 'lower');

P   = cholesky(pdMatrix);

check = roundn(pdMatrix, -10) == roundn(P * P', -10);
if sum(sum(check)) == n * n
    sprintf('Equal')
else
    sprintf('Not Equal')
end

Github → https://github.com/Hom-Wang/MATLAB/tree/master/cholesky

KitSprout

Cholesky 正定矩陣分解

Comments