矩阵求逆的简单记录

写在前面

话说第一次完整的编写矩阵求逆函数是在写结构力学求解器时候，那时候用的是初等变换法，效率较低，不过那时候还没顾得上效率，求解器没有开发完整，记得只支持连续跨简支梁求弯矩（不记得这个术语还准确不...），好吧，是时候引出正文了。

 

4X4转换矩阵求逆方法

求逆矩阵方法很多，效率不一，而研究矩阵数据结构的设计就可以是一个领域了，怎么样占用内存小？怎么样运算效率高，当然本文不涉及这些的研究，有兴趣可以自行搜索学习这方面知识。

https://zhuanlan.zhihu.com/p/50093662

https://www.cnblogs.com/johnson0522/archive/2009/06/09/1499515.html

常用的有待定系数法、初等变换法、伴随矩阵法，当然对于固定阶的矩阵可以直接推出逆矩阵行列元素的值。

先说下图形学中转换矩阵逆矩阵快速求法吧。注意此处的是矩阵乘列向量形式的转换矩阵求逆，即求逆的过程中利用了第三行是0 0 0 1的特点。

// 求4x4 矩阵的逆矩阵
// m 源矩阵   mi 源矩阵的逆矩阵
bool MatrixInverse(const Transform* trs, Transform* trsRe)
{
    const float* values = trs->values();

    // 矩阵的行列式
    float det = (values[0] * (values[4] * values[8] - values[7] * values[5]) -
        values[3] * (values[1] * values[8] - values[7] * values[2]) +
        values[6] * (values[1] * values[5] - values[4] * values[2]));

    // 先判断行列式是否为0。
    if (abs(det) < 1e-6)
        return false;

    float det_inv = 1.0f / det;

    array<float, 12> valuesRe;
    valuesRe[0] = det_inv * (values[4] * values[8] - values[7] * values[5]);
    valuesRe[3] = -det_inv * (values[3] * values[8] - values[6] * values[5]);
    valuesRe[6] = det_inv * (values[3] * values[7] - values[6] * values[4]);
    valuesRe[9] = -(valuesRe[0] * values[9] + valuesRe[3] * values[10] + valuesRe[6] * values[11]);

    valuesRe[1] = -det_inv * (values[1] * values[8] - values[7] * values[2]);
    valuesRe[4] = det_inv * (values[0] * values[8] - values[6] * values[2]);
    valuesRe[7] = -det_inv * (values[0] * values[7] - values[6] * values[1]);
    valuesRe[10] = -(valuesRe[1] * values[9] + valuesRe[4] * values[10] + valuesRe[7] * values[11]);

    valuesRe[2] = det_inv * (values[1] * values[5] - values[4] * values[2]);
    valuesRe[5] = -det_inv * (values[0] * values[5] - values[3] * values[2]);
    valuesRe[8] = det_inv * (values[0] * values[4] - values[3] * values[1]);
    valuesRe[11] = -(valuesRe[2] * values[9] + valuesRe[5] * values[10] + valuesRe[8] * values[11]);

    trsRe->setValues(&(valuesRe[0]));

    return true;
}

其实这里用的是3X3矩阵求逆的伴随矩阵法+待定系数法来求解4X4的转换矩阵逆矩阵。

伴随矩阵法：详细请打开线性代数书籍或搜索。

定理：n阶矩阵A=[a ij]

 

 

注意上述代码中矩阵的values是这样排列的，请对应换成自己的矩阵结构中的元素，

 初等变换法求NxN矩阵的逆矩阵

c++代码

void MatrixAXA::GetInverseMatrix(const MatrixAXA& matPara, MatrixAXA& matInverse)
{
    MatrixAXA mat(matPara);

    //  �����
    int nStep = mat.m_mat.size();
    matInverse = MatrixAXA(nStep);
    for (int i = 0; i < nStep; i++)
    {
        matInverse.m_mat[i][i] = 1.0;
    }

    //  �õ�һ�н��·��еĵ�һ�б�Ϊ0���ڶ��н��·��еĵڶ��б�Ϊ0...
    for (int i = 0; i < mat.m_mat.size() - 1; i++)
    {
        if (RealEq(mat.m_mat[i][i], 0.0))
        {
            //  �����i�еĵ�i��Ϊ0
            int nTemp = 0;
            for(auto ptrRow = mat.m_mat.begin(); ptrRow != mat.m_mat.end(); ++ptrRow)
            {
                if (!RealEq(mat.m_mat[ptrRow->first][i], 0.0) && ptrRow->first > i)
                {
                    nTemp = ptrRow->first;
                    break;
                }
            }

            //  ��i�и��м��ϸ��и��еõ��µĵ�i��
            for (int colum_i = 0; colum_i < mat.m_mat[i].size(); colum_i++)
            {
                mat.m_mat[i][colum_i] += mat.m_mat[nTemp][colum_i];               //  ******
                matInverse.m_mat[i][colum_i] += matInverse.m_mat[nTemp][colum_i]; //  ******
            }

        }

        //  �任Ŀ���е�i��Ϊ0
        for (int j = i + 1; j < mat.m_mat.size(); j++)
        {
            if (RealEq(mat.m_mat[j][i], 0.0) ||
                RealEq(mat.m_mat[i][i], 0.0))
                continue;

            double dTemp = -mat.m_mat[j][i] / mat.m_mat[i][i];
            for (int colum_i = 0; colum_i < mat.m_mat[j].size(); colum_i++)
            {
                mat.m_mat[j][colum_i] += mat.m_mat[i][colum_i] * dTemp;                //  ******
                matInverse.m_mat[j][colum_i] += matInverse.m_mat[i][colum_i] * dTemp;  //  ******
            }

        }
    }

    //  �����������һ�н�ǰ���еĵ�nStep�б任Ϊ0...
    for (int rowInv = nStep - 1; rowInv > 0; rowInv--)
    {
        for (int row_i = rowInv - 1; row_i >= 0; row_i--)
        {
            if (RealEq(mat.m_mat[row_i][rowInv], 0.0) ||
                RealEq(mat.m_mat[rowInv][rowInv], 0.0))
                continue;

            //  �任
            double dTemp = -mat.m_mat[row_i][rowInv] / mat.m_mat[rowInv][rowInv];
            for (int colum_i = 0; colum_i < mat.m_mat[row_i].size(); colum_i++)
            {
                mat.m_mat[row_i][colum_i] += dTemp * mat.m_mat[rowInv][colum_i];                //  ******
                matInverse.m_mat[row_i][colum_i] += dTemp * matInverse.m_mat[rowInv][colum_i];  //  ******
            }
        }

    }

    //  �����е�λ��
    for (int row_i = 0; row_i < nStep; row_i++)
    {
        double dTe1 = mat.m_mat[row_i][row_i];
        for (int colum_i = 0; colum_i < nStep; colum_i++)
        {
            if (RealEq(mat.m_mat[row_i][row_i], 0.0)) continue;
            mat.m_mat[row_i][colum_i] /= dTe1;                //  ******
            matInverse.m_mat[row_i][colum_i] /= dTe1;  //  ******
        }
    }  //  end for
}

 c#代码

/// <summary>
      /// 求逆矩阵
      /// 初等变换法
      /// </summary>
      public static void GetInverseMatrix(MatrixAXA matPara, ref MatrixAXA matInverse)
      {
         var mat = new MatrixAXA(matPara);

         //  单位矩阵
         int nStep = mat.m_mat.Count;
         matInverse = new MatrixAXA(nStep);
         for (int i = 0; i < nStep; i++)
         {
            matInverse.m_mat[i][i] = 1.0;
         }

         //  用第一行将下方行的第一列变为0，第二行将下方行的第二列变为0...
         for (int i = 0; i < mat.m_mat.Count - 1; i++)
         {
            if (MathTool.RealEq(mat.m_mat[i][i], 0.0))
            {
               //  如果第i行的第i列为0
               int nTemp = 0;
               foreach (var item in mat.m_mat.Keys)
               {
                  if (!MathTool.RealEq(mat.m_mat[item][i], 0.0) && item > i)
                  {
                     nTemp = item;
                     break;
                  }
               }

               //  第i行各列加上该行各列得到新的第i行
               for (int colum_i = 0; colum_i < mat.m_mat[i].Count; colum_i++)
               {
                  mat.m_mat[i][colum_i] += mat.m_mat[nTemp][colum_i];               //  ******
                  matInverse.m_mat[i][colum_i] += matInverse.m_mat[nTemp][colum_i]; //  ******
               }

            }

            //  变换目标行第i列为0
            for (int j = i + 1; j < mat.m_mat.Count; j++)
            {
               if (MathTool.RealEq(mat.m_mat[j][i], 0.0) ||
                 MathTool.RealEq(mat.m_mat[i][i], 0.0))
                  continue;

               double dTemp = -mat.m_mat[j][i] / mat.m_mat[i][i];
               for (int colum_i = 0; colum_i < mat.m_mat[j].Count; colum_i++)
               {
                  mat.m_mat[j][colum_i] += mat.m_mat[i][colum_i] * dTemp;                //  ******
                  matInverse.m_mat[j][colum_i] += matInverse.m_mat[i][colum_i] * dTemp;  //  ******
               }

            }
         }

         //  反过来，最后一行将前面行的第nStep列变换为0...
         for (int rowInv = nStep - 1; rowInv > 0; rowInv--)
         {
            for (int row_i = rowInv - 1; row_i >= 0; row_i--)
            {
               if (MathTool.RealEq(mat.m_mat[row_i][rowInv], 0.0) ||
                 MathTool.RealEq(mat.m_mat[rowInv][rowInv], 0.0))
                  continue;

               //  变换
               double dTemp = -mat.m_mat[row_i][rowInv] / mat.m_mat[rowInv][rowInv];
               for (int colum_i = 0; colum_i < mat.m_mat[row_i].Count; colum_i++)
               {
                  mat.m_mat[row_i][colum_i] += dTemp * mat.m_mat[rowInv][colum_i];                //  ******
                  matInverse.m_mat[row_i][colum_i] += dTemp * matInverse.m_mat[rowInv][colum_i];  //  ******
               }
            }

         }

         //  将各行单位化
         for (int row_i = 0; row_i < nStep; row_i++)
         {
            double dTe1 = mat.m_mat[row_i][row_i];
            for (int colum_i = 0; colum_i < nStep; colum_i++)
            {
               if (MathTool.RealEq(mat.m_mat[row_i][row_i], 0.0)) continue;
               mat.m_mat[row_i][colum_i] /= dTe1;                //  ******
               matInverse.m_mat[row_i][colum_i] /= dTe1;  //  ******
            }
         }  //  end for
      }

 C++和C#的都经过测试，不过代码是几年前写的，要参考的自行去重构吧~

总结

本文算是待定系数法、初等变换法、伴随矩阵法都有所介绍和应用，好了就到这里。