《pbrt》一書配套源碼上對于矩陣求逆使用的是全主元消去法，但是它的實現與我所見到過的全主元消去法還是略有不同的，有的地方還是值得思考一下的。

學過線性代數的應該都知道求矩陣的逆有伴隨矩陣法和初等變換法等方法，而高斯消元法、列主元消去法、全主元消去法這些算法都基于初等變換法，簡單說一下，對于矩陣$A$，我們想要求一個矩陣$B$，使得$AB=I$，那么可以先對$A$做一系列的初等行（列）變換，不妨做初等行變換，讓A變成單位矩陣，對于任何一個可逆矩陣來說，這是一定能夠做到的，即在$A$的左側乘上一系列的初等變換矩陣使得$P_1{P}_2{P}_3{P}_4???P_{i}A=I$，因此可以發現$P_1{P}_2{P}_3{P}_4???P_i=B=P_1{P}_2{P}_3{P}_4???P_{i}I$，那么我們可以在對A做初等行變換的同時對單位矩陣也做完全相同的初等行變換，最終得到的就應該是$A$的逆矩陣。

那么各種主元消去法就是基于這一點，全主元消去法的策略是每次選取矩陣中絕對值最大的一個元素為軸點，利用交換把它放到矩陣的最左上角，然后讓軸點所在列的其他元素利用初等行變換全部變為0，之后忽略掉軸點所在的行和列，此后的軸點選取都在前一個軸點的右下角的子矩陣內進行。

如果覺得還是不太明白的話可以去查閱各種消去法的具體算法。

這個是PBRT中的全主元消去法的實現：

Matrix4x4 Inverse(const Matrix4x4 &m) {int indxc[4], indxr[4];int ipiv[4] = {0, 0, 0, 0};Float minv[4][4];memcpy(minv, m.m, 4 * 4 * sizeof(Float));for (int i = 0; i < 4; i++) {int irow = 0, icol = 0;Float big = 0.f;// Choose pivotfor (int j = 0; j < 4; j++) {if (ipiv[j] != 1) {for (int k = 0; k < 4; k++) {if (ipiv[k] == 0) {if (std::abs(minv[j][k]) >= big) {big = Float(std::abs(minv[j][k]));irow = j;icol = k;}} else if (ipiv[k] > 1)Error("Singular matrix in MatrixInvert");}}}++ipiv[icol];// Swap rows _irow_ and _icol_ for pivotif (irow != icol) {for (int k = 0; k < 4; ++k) std::swap(minv[irow][k], minv[icol][k]);}indxr[i] = irow;indxc[i] = icol;if (minv[icol][icol] == 0.f) Error("Singular matrix in MatrixInvert");// Set $m[icol][icol]$ to one by scaling row _icol_ appropriatelyFloat pivinv = 1. / minv[icol][icol];minv[icol][icol] = 1.;for (int j = 0; j < 4; j++) minv[icol][j] *= pivinv;// Subtract this row from others to zero out their columnsfor (int j = 0; j < 4; j++) {if (j != icol) {Float save = minv[j][icol];minv[j][icol] = 0;for (int k = 0; k < 4; k++) minv[j][k] -= minv[icol][k] * save;}}}// Swap columns to reflect permutationfor (int j = 3; j >= 0; j--) {if (indxr[j] != indxc[j]) {for (int k = 0; k < 4; k++)std::swap(minv[k][indxr[j]], minv[k][indxc[j]]);}}return Matrix4x4(minv); }

關于這個實現，我覺得有以下幾個問題值得考慮：
1.這個實現是怎么做到不額外使用空間來存儲單位矩陣，然后對單位矩陣同時做初等變換得到逆矩陣，而是直接在原矩陣上做相應的變換得到逆矩陣的？
2.為什么在選定軸點后，假設軸點的位置為(i,j)，要交換第i行和第j行？
3.為什么在算法的結尾要以交換行時 逆過來的順序交換列？
4.這個實現是怎么判定奇異矩陣的？

讀者不妨可以自己先思考一下這幾個問題。

下面是我的一點想法：

1.整個算法確實完全都在對原矩陣進行操作，沒有任何對單位矩陣做的額外操作，但是其中有兩個地方值得注意：

首先是此處：

// Set $m[icol][icol]$ to one by scaling row _icol_ appropriatelyFloat pivinv = 1. / minv[icol][icol];minv[icol][icol] = 1.;for (int j = 0; j < 4; j++) minv[icol][j] *= pivinv;

可以看到，如果按照一般的全主元消去策略，這里是對軸點所在的那一行整體乘上一個系數，來讓軸點所在的位置變為1。但是此處多的這一行

minv[icol][icol] = 1.;

就是不對的，因為我們要讓原矩陣的主對角線上的元素變成1，但是加了這一句后，主對角線上的元素卻變成了pivinv，想一下如果我們對單位矩陣做相同的初等變換，主對角線上的元素確實應該是pivinv，因此這里實際上是讓原矩陣的主對角線的元素與單位矩陣經過相同初等變換后的對角線的元素相同，而非主對角線上的元素與原矩陣變換后的相同。

同樣的方式，在隨后讓軸點所在列的其他元素變為0時，也多了一行：

Float save = minv[j][icol];minv[j][icol] = 0;for (int k = 0; k < 4; k++) minv[j][k] -= minv[icol][k] * save; minv[j][icol] = 0;

因為如果對單位矩陣做這樣的變換，軸點所在列的其他元素均為0，因此應該是$0?=minv[icol][k]?save$，那么這里依舊是將單位矩陣變換后的數直接存放到原矩陣中，但是對于一行中的非軸點列的其他元素，不會受到影響。

這樣我們發現經過這樣的處理后軸點所在列的所有元素都與單位矩陣經過變換后相同列上的元素 完全相同，而非軸點所在列與原矩陣經過變換后的結果 完全相同，也就是說非軸點所在的列完全沒有受到這樣的處理的影響，不會對之后處理其他列做變換時使用的系數造成影響，并且同時保證了最終處理的結果就是單位矩陣經過相同變換后得到的結果。

2.在每次選定軸點(i,j)后，會將第i行和第j行互換，在全主元消去法中，我們每次都要讓軸點變換到當前處理矩陣的左上角，即主對角線上的位置，讓第i行和第j行互換，就使得軸點的位置變成了(j,j)，即第j個主對角元的位置，也就是說在這個實現中，并沒有每次讓軸點變到最左上角的主對角元處，而是讓它放到它所在的列的主對角元處。
然后在此后選取軸點時，就直接忽略掉第j列上的所有元素。

3.在算法結束時會讓矩陣的列按照行交換時逆著的順序來交換，這是因為算法結束時得到的矩陣應該是原矩陣的逆矩陣，而不是原矩陣經過一定的行交換后的逆矩陣。比如原來矩陣的第i行，與所得到的結果矩陣的第i列點乘的結果就應該是1，而不應該是與其他列點乘結果為1，如果在之前處理中我們讓第1行和第2行交換，又讓第2行與第4行交換，即$(1,2,3,4)?(2,1,3,4)?(2,4,3,1)$，那么最終得到的結果矩陣應該是：
第1列與原矩陣的第2行點乘結果為1
第2列與原矩陣的第4行點乘結果為1
第3列與原矩陣的第3行點乘結果為1
第4列與原矩陣的第1行點乘結果為1
這顯然是錯位了的，那么就需要對結果矩陣做列交換進行調整，并且應該以行交換的順序逆轉過來進行交換。

4.如果一個矩陣是奇異矩陣，在做初等行變換的過程中，一定會出現至少有一行為全0，這樣一定存在某個階段，軸點只能選0，沒有其他任何絕對值比0大的可供選擇的軸點，因此此時就可判定矩陣為奇異矩陣，不可逆。

總結

以上是生活随笔為你收集整理的pbrt源码中用全主元消去法求矩阵逆的实现的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。