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

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

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

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

這個(gè)是PBRT中的全主元消去法的實(shí)現(xiàn)：

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); }

關(guān)于這個(gè)實(shí)現(xiàn)，我覺得有以下幾個(gè)問(wèn)題值得考慮：
1.這個(gè)實(shí)現(xiàn)是怎么做到不額外使用空間來(lái)存儲(chǔ)單位矩陣，然后對(duì)單位矩陣同時(shí)做初等變換得到逆矩陣，而是直接在原矩陣上做相應(yīng)的變換得到逆矩陣的？
2.為什么在選定軸點(diǎn)后，假設(shè)軸點(diǎn)的位置為(i,j)，要交換第i行和第j行？
3.為什么在算法的結(jié)尾要以交換行時(shí) 逆過(guò)來(lái)的順序交換列？
4.這個(gè)實(shí)現(xiàn)是怎么判定奇異矩陣的？

讀者不妨可以自己先思考一下這幾個(gè)問(wèn)題。

下面是我的一點(diǎn)想法：

1.整個(gè)算法確實(shí)完全都在對(duì)原矩陣進(jìn)行操作，沒(méi)有任何對(duì)單位矩陣做的額外操作，但是其中有兩個(gè)地方值得注意：

首先是此處：

// 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;

可以看到，如果按照一般的全主元消去策略，這里是對(duì)軸點(diǎn)所在的那一行整體乘上一個(gè)系數(shù)，來(lái)讓軸點(diǎn)所在的位置變?yōu)?。但是此處多的這一行

minv[icol][icol] = 1.;

就是不對(duì)的，因?yàn)槲覀円屧仃嚨闹鲗?duì)角線上的元素變成1，但是加了這一句后，主對(duì)角線上的元素卻變成了pivinv，想一下如果我們對(duì)單位矩陣做相同的初等變換，主對(duì)角線上的元素確實(shí)應(yīng)該是pivinv，因此這里實(shí)際上是讓原矩陣的主對(duì)角線的元素與單位矩陣經(jīng)過(guò)相同初等變換后的對(duì)角線的元素相同，而非主對(duì)角線上的元素與原矩陣變換后的相同。

同樣的方式，在隨后讓軸點(diǎn)所在列的其他元素變?yōu)?時(shí)，也多了一行：

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;

因?yàn)槿绻麑?duì)單位矩陣做這樣的變換，軸點(diǎn)所在列的其他元素均為0，因此應(yīng)該是$0?=minv[icol][k]?save$，那么這里依舊是將單位矩陣變換后的數(shù)直接存放到原矩陣中，但是對(duì)于一行中的非軸點(diǎn)列的其他元素，不會(huì)受到影響。

這樣我們發(fā)現(xiàn)經(jīng)過(guò)這樣的處理后軸點(diǎn)所在列的所有元素都與單位矩陣經(jīng)過(guò)變換后相同列上的元素 完全相同，而非軸點(diǎn)所在列與原矩陣經(jīng)過(guò)變換后的結(jié)果 完全相同，也就是說(shuō)非軸點(diǎn)所在的列完全沒(méi)有受到這樣的處理的影響，不會(huì)對(duì)之后處理其他列做變換時(shí)使用的系數(shù)造成影響，并且同時(shí)保證了最終處理的結(jié)果就是單位矩陣經(jīng)過(guò)相同變換后得到的結(jié)果。

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

3.在算法結(jié)束時(shí)會(huì)讓矩陣的列按照行交換時(shí)逆著的順序來(lái)交換，這是因?yàn)樗惴ńY(jié)束時(shí)得到的矩陣應(yīng)該是原矩陣的逆矩陣，而不是原矩陣經(jīng)過(guò)一定的行交換后的逆矩陣。比如原來(lái)矩陣的第i行，與所得到的結(jié)果矩陣的第i列點(diǎn)乘的結(jié)果就應(yīng)該是1，而不應(yīng)該是與其他列點(diǎn)乘結(jié)果為1，如果在之前處理中我們讓第1行和第2行交換，又讓第2行與第4行交換，即$(1,2,3,4)?(2,1,3,4)?(2,4,3,1)$，那么最終得到的結(jié)果矩陣應(yīng)該是：
第1列與原矩陣的第2行點(diǎn)乘結(jié)果為1
第2列與原矩陣的第4行點(diǎn)乘結(jié)果為1
第3列與原矩陣的第3行點(diǎn)乘結(jié)果為1
第4列與原矩陣的第1行點(diǎn)乘結(jié)果為1
這顯然是錯(cuò)位了的，那么就需要對(duì)結(jié)果矩陣做列交換進(jìn)行調(diào)整，并且應(yīng)該以行交換的順序逆轉(zhuǎn)過(guò)來(lái)進(jìn)行交換。

4.如果一個(gè)矩陣是奇異矩陣，在做初等行變換的過(guò)程中，一定會(huì)出現(xiàn)至少有一行為全0，這樣一定存在某個(gè)階段，軸點(diǎn)只能選0，沒(méi)有其他任何絕對(duì)值比0大的可供選擇的軸點(diǎn)，因此此時(shí)就可判定矩陣為奇異矩陣，不可逆。

總結(jié)

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

如果覺得生活随笔網(wǎng)站內(nèi)容還不錯(cuò)，歡迎將生活随笔推薦給好友。