技術標籤：python 選取矩陣中某些行的數值

目前市面上的Julia教程很多，但大部分教程都是從0開始，對我這種有其他語言程式設計基礎的十分不友好，我真的不關注你的Unicode支援有多棒好嗎，類似這種東西在教程裡大講特講對我來說屬實沒必要。另一部分則是過於高深，各種高階特性滿天飛，本想著先拿來寫點東西，看到奇怪的語法當場被勸退。恰好最近在學數值分析，借這個機會剛好練習一下Julia的使用。我的使用僅關注功能的實現，在語法上僅限能跑通，對其他細節則不去關注。

本文的程式碼可以在這裡看到

先介紹一下本文用到的Julia語法

1. Julia的索引從1開始，與Python、C++不一樣
2. Array即我們需要的矩陣

x = [1. 2. 3.; 2 7 5; 1 4 9] # 3 * 3的矩陣
b = [1; 6; -3.] # 長度為3的列向量

3. size用來獲取矩陣的大小

size(x) # 返回二維矩陣的行數和列數
size(x, 1) # 返回矩陣的行數
size(x, 2) # 返回矩陣的列數

4. 矩陣的索引方法和matlab、python類似

x = [1 2 3; 4 5 6; 7 8 9]
x[1, 1] # 1
x[1, :] # 取第一行的值，返回一個列向量
x[end, end] # 9
x[:, 2] # 取第二列的值，返回一個列向量
 

5. 用.*可以實現對應元素相乘（和matlab一樣）

6. Julia的map和python類似，第一個引數為表示式，第二個引數為物件

高斯消去法

高斯消去法是用來求解線性方程組的經典解法，其基本思路是對增廣矩陣做行變換化為上三角矩陣後再回代求得方程組的解。以下均假設方程組有解。

首先定義函式的簽名

function gauss_elimination(X::Array, b::Array)

其中，X是n階方陣，b是n維列向量。

為了便於運算，我們將方陣和列向量合併為增廣矩陣

A = hcat(x, b)
row_num = size(A, 1)

接下來進行消去運算

for i in 1:row_num-1
    pivot = A[i, i] # 對角線上元素
    for j in i+1:row_num
        base = A[j, i] / pivot # 確定消去該行第一個非零元素的係數
        A[j,:] = A[j,:] - (base .* A[i,:]) # 對該行所有元素減去係數乘第i行的對應值
    end
end
 

以如下線性方程組為例

上述程式碼模擬的即是消去的過程

在得到變換後的增廣矩陣後即可進行回代

b = A[:, end]
b[end] /= A[end, end-1] # 最後一行直接計算即可
for i in row_num-1:-1:1 # 對其他行，先將除第一個非0列以外的其他列歸零
    pivot = A[i,i]
    b[i] -= sum(A[i,i+1:end-1] .* b[i+1:end])
    b[i] /= pivot
end

為了遵循課本的省記憶體的需求，將b作為結果返回。

高斯主元素法

高斯主元素法，又名列主元消去法，主要是為了解決在工程計算中由於計算機精度不足時除數過小所造成的精度丟失問題。其主要思想是，在每次消元時選取待選元素中最大的值作為當前目標。具體的思路可以在這裡看到。

主元素法的基本思路和普通的消去法一致，僅多了選擇最大列主元一項

target = map(item -> abs(item), A[i:end, i]) # 得到待選列
max_val = maximum(target) # 待選列最大值
max_index = findfirst(target .== max_val) # 待選列索引
if max_index != 1 # 若當前行的列值不是最大，則交換兩行
    A[i, :], A[max_index, :] = A[max_index, :], A[i, :]
end

參考文獻

1. 王開榮. 應用數值分析[M]. 高等教育出版社, 2010.