LDL分解并对厄而米特矩阵求逆（python仿真）

1. LDL分解

1.1. 原理

直接参考链接：https://blog.csdn.net/wanghu213/article/details/83963950这里不赘述。

1.2. python仿真

import numpy as np
np.set_printoptions(linewidth=np.inf)

M = 4

A = np.array([[3,      2+4j,   4+2.5j,   2+7.2j], 
              [2-4j,   4   ,   7+0.4j,   22+4j], 
              [4-2.5j, 7-0.4j, 5,        42+0.22j],
              [2-7.2j, 22-4j,  42-0.22j, 6]])

U = np.zeros_like(A)
L = np.zeros_like(A)
D = np.zeros_like(A)

for i in range(M):
    L[i, i] = 1

for j in range(M):
    tmp_sum = 0
    for m in range(j):
        tmp_sum += U[j, m] * L[j, m].conj()
    D[j, j] = A[j, j] - tmp_sum
    tmp_sum = 0
    for i in range(j+1, M):
        tmp_sum = 0
        for m in range(j):
            tmp_sum += U[i, m] * L[j, m].conj()
        U[i, j] = A[i, j] - tmp_sum
        L[i, j] = U[i, j] / D[j, j]

print('L:')
print(L)

print('D:')
print(D)

print('origin:')
print(A)

print('LDLH:')
print(np.dot(L, np.dot(D, L.conj().T)))

结果：

L:
[[ 1.        +0.j          0.        +0.j          0.        +0.j          0.        +0.j        ]
 [ 0.66666667-1.33333333j  1.        +0.j          0.        +0.j          0.        +0.j        ]
 [ 1.33333333-0.83333333j -0.375     +1.525j       1.        +0.j          0.        +0.j        ]
 [ 0.66666667-2.4j        -4.15      +0.7j         9.69471154+5.74278846j  1.        +0.j        ]]
D:
[[   3.        +0.00000000e+00j    0.        +0.00000000e+00j    0.        +0.00000000e+00j    0.        +0.00000000e+00j]
 [   0.        +0.00000000e+00j   -2.66666667+0.00000000e+00j    0.        +0.00000000e+00j    0.        +0.00000000e+00j]
 [   0.        +0.00000000e+00j    0.        +0.00000000e+00j    4.16      +0.00000000e+00j    0.        +0.00000000e+00j]
 [   0.        +0.00000000e+00j    0.        +0.00000000e+00j    0.        +0.00000000e+00j -493.56293269+2.84217094e-14j]]
origin:
[[ 3.+0.j    2.+4.j    4.+2.5j   2.+7.2j ]
 [ 2.-4.j    4.+0.j    7.+0.4j  22.+4.j  ]
 [ 4.-2.5j   7.-0.4j   5.+0.j   42.+0.22j]
 [ 2.-7.2j  22.-4.j   42.-0.22j  6.+0.j  ]]
LDLH:
[[ 3.+0.00000000e+00j  2.+4.00000000e+00j  4.+2.50000000e+00j  2.+7.20000000e+00j]
 [ 2.-4.00000000e+00j  4.+0.00000000e+00j  7.+4.00000000e-01j 22.+4.00000000e+00j]
 [ 4.-2.50000000e+00j  7.-4.00000000e-01j  5.+0.00000000e+00j 42.+2.20000000e-01j]
 [ 2.-7.20000000e+00j 22.-4.00000000e+00j 42.-2.20000000e-01j  6.+2.84217094e-14j]]**加粗样式**

可以看出，原始矩阵和$LD{L}^H$结果一致。

2. Hermitian矩阵求逆

主要思路是已知：
$LD{L}^H{A}^{-1}=I$; 设$D{L}^H{A}^{-1}=B$，则$LB=I$; 求解$B$， 其

2.1. 仿真代码

import numpy as np
import time
np.set_printoptions(linewidth=np.inf)
np.set_printoptions(precision=8,suppress=True)

M = 4

def LDLH(A):
    U = np.zeros_like(A)
    L = np.zeros_like(A)
    D = np.zeros_like(A)

    for i in range(M):
        L[i, i] = 1

    for j in range(M):
        tmp_sum = 0
        for m in range(j):
            tmp_sum += U[j, m] * L[j, m].conj()
        D[j, j] = A[j, j] - tmp_sum
        tmp_sum = 0
        for i in range(j+1, M):
            tmp_sum = 0
            for m in range(j):
                tmp_sum += U[i, m] * L[j, m].conj()
            U[i, j] = A[i, j] - tmp_sum
            L[i, j] = U[i, j] / D[j, j]

    return L, D

def LDL_solve(L, D):
    M = L.shape[0]
    B = np.zeros_like(L)
    V = np.zeros_like(L)
    for i in range(M):
        B[i, i] = 1
    for j in range(M-1):
        B[j+1, j] = -L[j+1, j]
        for i in range(j+2, M):
            tmp_sum = 0
            for m in range(j+1, i):
                tmp_sum += B[m, j] * L[i, m]
            B[i, j] = -L[i, j] - tmp_sum

    for j in range(M):
        V[M-1, j] = B[M-1, j] / D[M-1, M-1]
        V[j, M-1] = V[M-1, j].conj()
        for i in range(M-2, j-1, -1):
            tmp_sum = 0
            for m in range(i+1, M):
                tmp_sum += L[m, i].conj() * V[m, j]
            V[i, j] = B[i, j] / D[i, i] - tmp_sum
            V[j, i] = V[i, j].conj()

    return V

A = np.array([[3,      2+4j,   4+2.5j,   2+7.2j], 
              [2-4j,   4   ,   7+0.4j,   22+4j], 
              [4-2.5j, 7-0.4j, 5,        42+0.22j],
              [2-7.2j, 22-4j,  42-0.22j, 6]])


V_np = np.linalg.inv(A)

L, D = LDLH(A)
V = LDL_solve(L, D)


print('v_np - v:')
print(V_np - V)

print('AV')
print(np.dot(A, V))

结果：

v_np - v:
[[-0.+0.j  0.-0.j -0.+0.j  0.+0.j]
 [ 0.+0.j -0.-0.j -0.-0.j -0.+0.j]
 [ 0.-0.j -0.+0.j -0.+0.j  0.-0.j]
 [-0.-0.j -0.-0.j  0.+0.j -0.+0.j]]
AV
[[ 1.-0.j  0.+0.j -0.+0.j  0.-0.j]
 [-0.-0.j  1.-0.j  0.+0.j -0.-0.j]
 [-0.-0.j  0.-0.j  1.+0.j -0.-0.j]
 [-0.+0.j  0.-0.j  0.+0.j  1.-0.j]]

可以看出求逆的结果和numpy的计算结果一致，且原始矩阵和求逆矩阵乘积为单位阵。