计算 n x n 矩阵的倒数

Calculate the inverse of an n x n matrix

提问人:陳逸軒 提问时间:6/8/2023 最后编辑:陳逸軒 更新时间:6/9/2023 访问量:192

问:

我想使用找到佐剂矩阵和行列式的方法找到 n 阶矩阵的逆矩阵,然后推导出逆矩阵。我使用的是 C++ 语言,为了计算行列式,我使用了递归的概念。我想知道是否有更好的方法。

p.s. 由于分配的限制,不允许使用其他库。

矩阵.h

class Matrix
{
private:
  friend ostream &operator<<(ostream &, const Matrix &);
  friend istream &operator>>(istream &, Matrix &);

  friend Matrix getADJ(const Matrix);
  friend float getDet(const Matrix);

public:
  Matrix();
  Matrix(const int, const int, const float = 0.0f);
  Matrix(const Matrix &);
  ~Matrix();

  // operation
  Matrix operator/(const float) const;
  Matrix &operator=(const Matrix &);
  Matrix inverse() const;

private:
  float **data;
  int row;
  int column;
};

矩阵.cpp

#include "Matrix.h"

// self define func.
static float getRandFloat(const float lower = 0.0F, const float upper = 10.0F)
{
  float temp = (float)rand() / RAND_MAX;
  return lower + temp * (upper - lower);
}
static float **new_2DArray(const int _r, const int _c, const float init_val = 0.0f)
{
  if (_r < 0 || _c < 0)
  {
    cout << "[ERROR] Size should be a positive number!!" << endl;

    return nullptr;
  }
  const unsigned long r_ul = (unsigned long)_r;
  const unsigned long c_ul = (unsigned long)_c;

  float **temp = new float *[r_ul];
  for (unsigned long r = 0; r < r_ul; r++)
    temp[r] = new float[c_ul];

  for (unsigned long r = 0; r < r_ul; r++)
    for (unsigned long c = 0; c < c_ul; c++)
      temp[r][c] = init_val;

  return temp;
}

// matrix constructor
Matrix::Matrix()
{
  data = new_2DArray(row = 2, column = 2);
  // set row to 2, column to 2
  // all values in matrix are set to default value 0.0f
}
Matrix::Matrix(const int row, const int col, const float init_val)
{
  data = new_2DArray(this->row = row, this->column = col, init_val);
}
Matrix::Matrix(const Matrix &obj)
{
  row = obj.row;
  column = obj.column;
  data = new_2DArray(row, column);
  // the new matrix has the same value for row and column
  // create a new data space for the new matrix, assign all value equals to matrix obj

  for (int r = 0; r < row; r++)
    for (int c = 0; c < column; c++)
      data[r][c] = obj.data[r][c];
}
Matrix::~Matrix()
{
  for (int i = 0; i < row; i++)
    delete[] data[i];

  delete data;
}

//operators
Matrix &Matrix::operator=(const Matrix &obj)
{
  if (row != obj.row || column != obj.column)
  {
    cout << "[ERROR] Invalid operation: Size different!!" << endl;
    return *this;
  }
  Matrix *temp = new Matrix(obj);

  data = temp->data;

  return *this;
}

Matrix Matrix::inverse() const
{
  if (row != column || row <= 1)
  {
    cout << "[ERROR] Invalid operation: should be a nxn matrix(n >= 2)!!" << endl;
    return *this;
  }
  Matrix temp(*this);
  float det = getDet(temp);
  if (det == 0.0f)
  {
    cout << "[ERROR] DET should`t be zero!!" << endl;
    return *this;
  }

  return getADJ(temp) / det;
};

Matrix getADJ(const Matrix m)
{
  const int n = m.row;
  Matrix result(n, n);

  for (int c = 0; c < n; c++)
  {
    const int limited_column = c;
    bool if_positive = (c % 2 ? false : true);
    for (int r = 0; r < n; r++)
    {
      const int limited_row = r;
      Matrix temp(n - 1, n - 1);

      for (int sr = 0; sr < n; sr++)
      {
        if (sr == limited_row)
          continue;
        for (int sc = 0; sc < n; sc++)
        {
          if (sc == limited_column)
            continue;

          temp.data[sr > limited_row ? sr - 1 : sr][sc > limited_column ? sc - 1 : sc] = m.data[sr][sc];
        }
      }
      result.data[c][r] = getDet(temp) * (if_positive ? 1.0f : -1.0f);
      if_positive = !if_positive;
    }
  }

  return result;
}

float getDet(const Matrix m)
{
  const int n = m.row;
  if (n == 1)
    return m.data[0][0];

  float curLevel = 0.0f;
  for (int i = 0; i < n; i++)
  {
    const int limited_column = i;
    // limit row is always equal to zero, when calculating DET
    Matrix temp(n - 1, n - 1);
    for (int r = 1; r < n; r++)
    {
      for (int c = 0; c < n; c++)
      {
        if (c == limited_column)
          continue;
        temp.data[r - 1][c > limited_column ? c - 1 : c] = m.data[r][c];
      }

      curLevel += (m.data[0][i] * getDet(temp) * (i % 2 ? -1.0f : 1.0f));
    }
  }

  return curLevel;
}

// friend function: overloading << and >>
ostream &operator<<(ostream &out, const Matrix &m)
{
  for (int r = 0; r < m.row; r++)
  {
    for (int c = 0; c < m.column; c++)
      out << setw(10) << fixed << setprecision(5) << m.data[r][c];
    if (r != m.row - 1)
      out << endl;
  }

  return out;
}
istream &operator>>(istream &in, Matrix &m)
{
  for (int r = 0; r < m.row; r++)
    for (int c = 0; c < m.column; c++)
      m.data[r][c] = getRandFloat();
  return in;
}

main.cpp

#include <iostream>
#include "Matrix.h"

int main()
{
  srand((unsigned)time(0)); // set random seed

  Matrix m(6, 6); //m is a 6x6 matrix with all value set to 0
  std::cin >> m; // give each element in the matrix a random float value, when using >> operator

  cout << "inverse of m:" << endl;
  cout << m.inverse() << endl;
  //"<<" operator is use to output all value int the matrix
 
  return 0;
}

我目前正在测试的方法似乎工作正常。但是,我想确认它是否是最有效的解决方案。

输入

8.0867 3.15236 1.66545 1.26714 6.88726 4.24564
6.41356 2.64673 3.65965 7.65898 4.52873 4.3264
3.72865 7.36083 3.49968 9.13076 0.72390 6.62454
8.59957 2.96699 6.25450 9.29979 1.59475 2.93691
0.57955 0.50900 4.75473 2.79461 8.99372 7.43932
2.59571 6.15638 0.20169 9.80159 5.25058 6.5434

输出

0.01317 0.64093 0.12824 -0.34249 -0.13077 -0.25974
0.24609 -1.50393 -0.20641 0.81144 0.13351 0.52768
0.07594 -0.94817 -0.13510 0.61556 0.19399 0.21758
-0.10456 0.13971 -0.06629 0.002489 -0.02853 0.07390
0.13616 -0.92044 -0.32327 0.55861 0.14560 0.43126
-0.19173 1.71926 0.50620 -1.09852 -0.15382 -0.70406
C++ 递归 数学 -逆矩阵

评论

1赞 Pepijn Kramer 6/8/2023
不,您仍然执行过多的动态内存分配。使矩阵成为具有固定大小的模板类。话虽如此,如果你真的想编写自己的矩阵类,你唯一应该做的就是通过运行大量矩阵操作并从那里找到瓶颈来分析你的代码(它们可以是特定于平台的:依赖于缓存、分支预测、编译器的矢量化等)。总而言之,如果你真的需要快速的矩阵运算,请使用库
1赞 Pepijn Kramer 6/8/2023
请注意,递归通常也不是生成最快代码的最佳方法,更改为迭代算法(是的,每个递归算法都可以转换为迭代算法,或者相反)。
1赞 Pepijn Kramer 6/8/2023
作业,我不能使用外部库“部分是您应该在问题中提到的一大限制。另一个优化点,如果可以的话,避免循环中的 if 语句(多做一点计算甚至可以加快速度),错过分支预测也会减慢你的执行速度。
1赞 PaulMcKenzie 6/8/2023
您的类违反了 3 规则,因为它没有用户定义的赋值运算符。因此,任何与 的复制语义有关的事情,例如按值传递和返回,都是可疑的。 -- 如果析构函数尝试将内存 用于 ,则现在会出现双重删除错误。MatrixMatrixMatrix{Matrix m(6, 6); Matrix m2(5,5); m = m2;}delete[]Matrix
1赞 Spektre 6/8/2023
行列式基于逆的逆矩阵仅对小矩阵快速,对于较大的矩阵,使用 GEM(高斯消除法),它要快得多,但是由于精度问题,要使其鲁棒有点棘手......但是,它可以在没有任何递归的情况下就地完成,请参阅高斯消元法,没有加速结果,以获得一些灵感,GEM 和次行列式版本都有,甚至模板版本

答: 暂无答案