pbrt源码中用全主元消去法求矩阵逆的实现

《pbrt》一书配套源码上对于矩阵求逆使用的是全主元消去法，但是它的实现与我所见到过的全主元消去法还是略有不同的，有的地方还是值得思考一下的。

学过线性代数的应该都知道求矩阵的逆有伴随矩阵法和初等变换法等方法，而高斯消元法、列主元消去法、全主元消去法这些算法都基于初等变换法，简单说一下，对于矩阵AAA，我们想要求一个矩阵BBB，使得AB=IAB=IAB=I，那么可以先对AAA做一系列的初等行（列）变换，不妨做初等行变换，让A变成单位矩阵，对于任何一个可逆矩阵来说，这是一定能够做到的，即在AAA的左侧乘上一系列的初等变换矩阵使得P1P2P3P4⋅⋅⋅PiA=IP_1P_2P_3P_4···P_iA=IP1P2P3P4⋅⋅⋅PiA=I，因此可以发现P1P2P3P4⋅⋅⋅Pi=B=P1P2P3P4⋅⋅⋅PiIP_1P_2P_3P_4···P_i=B=P_1P_2P_3P_4···P_iIP1P2P3P4⋅⋅⋅Pi=B=P1P2P3P4⋅⋅⋅PiI，那么我们可以在对A做初等行变换的同时对单位矩阵也做完全相同的初等行变换，最终得到的就应该是AAA的逆矩阵。

那么各种主元消去法就是基于这一点，全主元消去法的策略是每次选取矩阵中绝对值最大的一个元素为轴点，利用交换把它放到矩阵的最左上角，然后让轴点所在列的其他元素利用初等行变换全部变为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]∗save0-=minv[icol][k] * save0−=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,3,4)-(2,1,3,4)-(2,4,3,1)(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源码中用全主元消去法求矩阵逆的实现相关推荐

NeuCF源码中用到的模块（函数）
论文:<Neural Collaborative Filtering>源码中用到的模块(函数) from keras.layers import Embedding, Input, Den ...
抖音seo源码，抖音seo矩阵系统源码技术搭建
抖音seo源码,抖音seo矩阵系统源码技术搭建抖音seo源码,抖音seo矩阵系统底层框架上支持了从ai视频混剪,视频批量原创产出,云端数字人视频制作,账号矩阵,视频一键分发,站内实现搜索排名,到同城 ...
抖音seo源码，抖音seo矩阵系统源码搭建技术+二开开源代码
抖音seo源码,抖音seo矩阵系统底层框架上支持了从ai视频创意混剪,上传素材批量产出,5个定时投放入口,账号矩阵搭建打造,视频内容优化,搜索引擎实现搜索排名,到同城客户裂变营销智能企业号营销管理4个 ...
抖音seo源码搭建，抖音矩阵系统源码分发，抖音矩阵同步分发
前言:抖音seo源码,抖音矩阵系统源码搭建,抖音矩阵同步分发.抖音seo源码部署是需要对接到这些正规接口再来做开发的,目前账号矩阵程序开发的功能,围绕一键管理多个账号,做到定时投放,关键词自动化生成霸 ...
抖音seo源码，抖音矩阵系统源码搭建，抖音矩阵同步分发。
前言:抖音seo源码,抖音矩阵系统源码搭建,抖音矩阵同步分发.抖音seo源码部署是需要对接到这些正规接口再来做开发的,目前账号矩阵程序开发的功能,围绕一键管理多个账号,做到定时投放,关键词自动化生成霸 ...
抖音seo源码，抖音seo矩阵系统分发源码技术搭建
抖音seo源码,抖音seo矩阵系统源码技术搭建,抖音seo源码技术开发思路梳理搭建抖音seo源码,抖音seo矩阵系统底层框架上支持了从ai视频混剪,视频批量原创产出,云存储批量视频制作,账号矩阵,视 ...
抖音SEO，抖音seo源码，抖音seo矩阵，抖音seo优化
抖音SEO,抖音seo源码,抖音seo矩阵,抖音seo优化抖音seo源码,抖音seo矩阵系统底层框架上支持了从ai视频混剪,视频批量原创产出,云端数字人视频制作,账号矩阵,视频一键分发,站内实现搜索 ...
抖音SEO，抖音seo源码，抖音seo矩阵，抖音seo独立部署
抖音SEO,抖音seo源码,抖音seo矩阵,抖音seo独立部署抖音SEO矩阵系统源码如何去做独立部署,首先我们要深入了解到这套系统的开发逻辑是什么呢?开发的前序是现在是抖音平台的流量激增,现在抖音站 ...
抖音矩阵系统源码搭建，抖音矩阵系统开发逻辑原理，抖音矩阵系统如何搭建？
抖音矩阵系统源码搭建,抖音矩阵系统开发逻辑原理,抖音矩阵系统如何搭建? 首先这里我们普及一下抖音矩阵系统大致可以分为四种类型,分别是:内容创作者矩阵.MCN机构/内容服务商矩阵.企业主/品牌主矩阵.媒 ...
抖音账号矩阵系统/抖音seo霸屏系统源码/关键词短视频账号矩阵源码/独立私有部署/可定制开发
前言:抖音账号矩阵系统/抖音seo霸屏系统源码/关键词短视频账号矩阵源码/独立私有部署/可定制开发场景:抖音账号矩阵系统/抖音seo霸屏系统/抖音矩阵seo系统源码/独立部署,技术团队如何围绕抖音矩 ...

pbrt源码中用全主元消去法求矩阵逆的实现

pbrt源码中用全主元消去法求矩阵逆的实现相关推荐

最新文章

热门文章