列主元方法是一种用于求解矩阵逆的数值方法，特别适用于在计算机上实现。其基本思想是通过高斯消元法将矩阵转换为上三角矩阵，然后通过回代求解矩阵的逆。以下是列主元方法求解矩阵 A A A 的逆的步骤：
步骤 1：初始化
构造增广矩阵
[
A
∣
I
]
[A | I]
[A∣I]，其中
I
I
I 是
n
n
n 阶单位矩阵。
步骤 2：列主元选择
对于第
k
k
k 列（
k
=
1
,
2
,
…
,
n
k = 1, 2, \ldots, n
k=1,2,…,n），找到列主元，即找到
i
k
i_k
ik 使得：
∣
a
i
k
,
k
∣
=
max
i
≥
k
∣
a
i
,
k
∣
|a_{i_k,k}| = \max_{i \geq k} |a_{i,k}|
∣aik,k∣=i≥kmax∣ai,k∣
如果
i
k
≠
k
i_k \neq k
ik=k，则交换第
k
k
k 行和第
i
k
i_k
ik 行。
步骤 3：高斯消元
对于每一列
k
k
k（
k
=
1
,
2
,
…
,
n
−
1
k = 1, 2, \ldots, n-1
k=1,2,…,n−1），进行以下操作：
- 归一化第
k
k
k 行的列主元：
a k , k ← 1 a k , k a_{k,k} \leftarrow \frac{1}{a_{k,k}} ak,k←ak,k1 - 更新第
k
k
k 行的其他元素：
a k , j ← a k , j a k , k 对于所有 j ≠ k a_{k,j} \leftarrow \frac{a_{k,j}}{a_{k,k}} \quad \text{对于所有 } j \neq k ak,j←ak,kak,j对于所有 j=k - 消去下方所有行的第
k
k
k 列元素：
对于所有 i > k i > k i>k，计算：
m i , k = a i , k m_{i,k} = a_{i,k} mi,k=ai,k
然后更新第 i i i 行：
a i , j ← a i , j − m i , k ⋅ a k , j 对于所有 j a_{i,j} \leftarrow a_{i,j} - m_{i,k} \cdot a_{k,j} \quad \text{对于所有 } j ai,j←ai,j−mi,k⋅ak,j对于所有 j
步骤 4：回代求解
当矩阵 A A A 被转换为上三角矩阵后，从最后一行开始回代：
对于每一行 k k k（ k = n , n − 1 , … , 1 k = n, n-1, \ldots, 1 k=n,n−1,…,1），进行以下操作： - 归一化第
k
k
k 行的最后一个非零元素（即对角线元素）：
a k , k ← 1 a k , k a_{k,k} \leftarrow \frac{1}{a_{k,k}} ak,k←ak,k1 - 更新第
k
k
k 行的其他元素：
a k , j ← a k , j a k , k 对于所有 j ≠ k a_{k,j} \leftarrow \frac{a_{k,j}}{a_{k,k}} \quad \text{对于所有 } j \neq k ak,j←ak,kak,j对于所有 j=k - 消去上方所有行的第
k
k
k 列元素：
对于所有 i < k i < k i<k，计算：
m i , k = a i , k m_{i,k} = a_{i,k} mi,k=ai,k
然后更新第 i i i 行：
a i , j ← a i , j − m i , k ⋅ a k , j 对于所有 j a_{i,j} \leftarrow a_{i,j} - m_{i,k} \cdot a_{k,j} \quad \text{对于所有 } j ai,j←ai,j−mi,k⋅ak,j对于所有 j
步骤 5：结果提取
经过上述步骤后，增广矩阵的左侧变为单位矩阵，而右侧则变为了 A A A 的逆矩阵 A − 1 A^{-1} A−1。提取右侧的矩阵即为所求的逆矩阵。
需要注意的是，上述公式中的 a i , j a_{i,j} ai,j 表示增广矩阵中的元素，包括原矩阵 A A A 和单位矩阵 I I I 的部分。在实际计算中，这些操作会同时应用于 A A A 和 I I I，最终 I I I 的位置会被 A − 1 A^{-1} A−1 所取代。
此外，如果在任何步骤中发现对角线上的元素 a k , k a_{k,k} ak,k 为零或非常接近零，那么矩阵 A A A 可能是奇异矩阵，无法求逆。在这种情况下，算法应该停止并报错。
Julia 代码
美化数据格式
using DataFrames
function pm(A,b)
m,n=size(A); z=[]
for i=1:n
st=i
z=[z; "a:$st"]
end
for i=1:n
st=i
z=[z;"b:$st"]
end
println(DataFrame([A b],z))
end
function luInv(A,par=false)
n=size(A,1);T=typeof(A[1,1])
A=copy(A); E = zeros(T,n,n);
for i=1:n E[i,i]=1//1 end
if par pm(A, E) end
if par println("化为上三角") end
for i=1:n-1
id=argmax(abs.(A[i:n,i])) # 寻找列主元
id=id-1
A[i,i:n], A[i+id,i:n]= A[i+id,i:n],A[i,i:n]
E[i,:], E[i+id,:] =E[i+id,:], E[i,:]
for j=i+1:n
c=A[j,i]/A[i,i]
E[j,:]=E[j,:]-E[i,:]*c
A[j,i:n]=A[j,i:n]-A[i,i:n]*c
end
if par pm(A, E) end
end
if par println("化为对角") end
for i=n:-1:2
for j=1:i-1
c=A[j,i]/A[i,i]
E[j,:]=E[j,:]-E[i,:]*c
A[j,1:i]=A[j,1:i]-A[i,1:i]*c
end
if par pm(A, E) end
end
IA=copy(E);
for j=1:n
IA[j,:] = E[j,:]/A[j,j]; A[j,j]=A[j,j]/A[j,j]
end
if par pm(A, IA) end
return(IA)
end
举例
n=3;
A=zeros(Rational,n,n)
for i=1:n-1
A[i,i]=0;
A[i,i+1]=11//1;
A[i+1,i]=7//1;
end
A[n,n]=3//1;
IA=luInv(A,true)
结果
