列主元方法是一種用于求解矩陣逆的數值方法，特別適用于在計算機上實現。其基本思想是通過高斯消元法將矩陣轉換為上三角矩陣，然后通過回代求解矩陣的逆。以下是列主元方法求解矩陣 $A$ 的逆的步驟：

步驟 1：初始化
構造增廣矩陣 $[A|I]$，其中 $I$ 是 $n$ 階單位矩陣。
步驟 2：列主元選擇
對于第 $k$ 列（$k=1,2,\dots ,n$），找到列主元，即找到 $i_k$ 使得：
$$|a_{i_k,k}|=\underset{i\ge k}{\max }|a_{i,k}|$$
如果 $i_k\ne k$，則交換第 $k$ 行和第 $i_k$ 行。
步驟 3：高斯消元
對于每一列 $k$（$k=1,2,\dots ,n?1$），進行以下操作：

• 歸一化第 $k$ 行的列主元：
  $$a_{k,k}\leftarrow \frac{1}{a_{k,k}}$$
• 更新第 $k$ 行的其他元素：
  $$a_{k,j}\leftarrow \frac{a_{k,j}}{a_{k,k}}\text{對于所有?}j\ne k$$
• 消去下方所有行的第 $k$ 列元素：
  對于所有 $i>k$，計算：
  $$m_{i,k}=a_{i,k}$$
  然后更新第 $i$ 行：
  $$a_{i,j}\leftarrow a_{i,j}?m_{i,k}?a_{k,j}\text{對于所有?}j$$
  步驟 4：回代求解
  當矩陣 $A$ 被轉換為上三角矩陣后，從最后一行開始回代：
  對于每一行 $k$（$k=n,n?1,\dots ,1$），進行以下操作：
• 歸一化第 $k$ 行的最后一個非零元素（即對角線元素）：
  $$a_{k,k}\leftarrow \frac{1}{a_{k,k}}$$
• 更新第 $k$ 行的其他元素：
  $$a_{k,j}\leftarrow \frac{a_{k,j}}{a_{k,k}}\text{對于所有?}j\ne k$$
• 消去上方所有行的第 $k$ 列元素：
  對于所有 $i<k$，計算：
  $$m_{i,k}=a_{i,k}$$
  然后更新第 $i$ 行：
  $$a_{i,j}\leftarrow a_{i,j}?m_{i,k}?a_{k,j}\text{對于所有?}j$$
  步驟 5：結果提取
  經過上述步驟后，增廣矩陣的左側變為單位矩陣，而右側則變為了 $A$ 的逆矩陣 $A^{?1}$。提取右側的矩陣即為所求的逆矩陣。
  需要注意的是，上述公式中的 $a_{i,j}$ 表示增廣矩陣中的元素，包括原矩陣 $A$ 和單位矩陣 $I$ 的部分。在實際計算中，這些操作會同時應用于 $A$ 和 $I$，最終 $I$ 的位置會被 $A^{?1}$ 所取代。
  此外，如果在任何步驟中發現對角線上的元素 $a_{k,k}$ 為零或非常接近零，那么矩陣 $A$ 可能是奇異矩陣，無法求逆。在這種情況下，算法應該停止并報錯。

Julia 代碼

美化數據格式

using DataFrames
function pm(A,b)m,n=size(A); z=[]for i=1:n  st=iz=[z; "a:$st"]endfor i=1:nst=iz=[z;"b:$st"]endprintln(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  endif par pm(A, E) endif par println("化為上三角")  endfor i=1:n-1id=argmax(abs.(A[i:n,i])) # 尋找列主元 id=id-1A[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:nc=A[j,i]/A[i,i]E[j,:]=E[j,:]-E[i,:]*cA[j,i:n]=A[j,i:n]-A[i,i:n]*cendif par pm(A, E) endendif par println("化為對角") endfor i=n:-1:2for j=1:i-1c=A[j,i]/A[i,i]E[j,:]=E[j,:]-E[i,:]*cA[j,1:i]=A[j,1:i]-A[i,1:i]*cendif par pm(A, E) endendIA=copy(E);for j=1:nIA[j,:] = E[j,:]/A[j,j]; A[j,j]=A[j,j]/A[j,j]endif par pm(A, IA) endreturn(IA)
end

舉例

n=3;
A=zeros(Rational,n,n)
for i=1:n-1A[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)

結果