简单的 3x3 矩阵逆代码 (C++)
声明：本页面是StackOverFlow热门问题的中英对照翻译，遵循CC BY-SA 4.0协议，如果您需要使用它，必须同样遵循CC BY-SA许可，注明原文地址和作者信息，同时你必须将它归于原作者（不是我）：StackOverFlow
原文地址: http://stackoverflow.com/questions/983999/
Warning: these are provided under cc-by-sa 4.0 license. You are free to use/share it, But you must attribute it to the original authors (not me):
StackOverFlow
Simple 3x3 matrix inverse code (C++)
提问by batty
What's the easiest way to compute a 3x3 matrix inverse?
计算 3x3 矩阵逆的最简单方法是什么？
I'm just looking for a short code snippet that'll do the trick for non-singular matrices, possibly using Cramer's rule. It doesn't need to be highly optimized. I'd prefer simplicity over speed. I'd rather not link in additional libraries.
我只是在寻找一个简短的代码片段，它可以解决非奇异矩阵的问题，可能使用 Cramer 规则。它不需要高度优化。我更喜欢简单而不是速度。我宁愿不链接其他库。
采纳答案by Suvesh Pratapa
Why don't you try to code it yourself? Take it as a challenge. :)
你为什么不尝试自己编码？把它当作一个挑战。:)
For a 3×3 matrix
对于 3×3 矩阵
(source: wolfram.com)
（来源：wolfram.com）
the matrix inverse is
矩阵逆是
(source: wolfram.com)
（来源：wolfram.com）
I'm assuming you know what the determinant of a matrix |A| is.
我假设你知道矩阵的行列式 |A| 是。
Images (c) Wolfram|Alphaand mathworld.wolfram(06-11-09, 22.06)
图片 (c) Wolfram|Alpha和 mathworld.wolfram(06-11-09, 22.06)
回答by Cornstalks
Here's a version of batty's answer, but this computes the correctinverse. batty's version computes the transpose of the inverse.
这是 batty 的答案的一个版本，但这计算了正确的逆。batty 的版本计算逆的转置。
// computes the inverse of a matrix m
double det = m(0, 0) * (m(1, 1) * m(2, 2) - m(2, 1) * m(1, 2)) -
m(0, 1) * (m(1, 0) * m(2, 2) - m(1, 2) * m(2, 0)) +
m(0, 2) * (m(1, 0) * m(2, 1) - m(1, 1) * m(2, 0));
double invdet = 1 / det;
Matrix33d minv; // inverse of matrix m
minv(0, 0) = (m(1, 1) * m(2, 2) - m(2, 1) * m(1, 2)) * invdet;
minv(0, 1) = (m(0, 2) * m(2, 1) - m(0, 1) * m(2, 2)) * invdet;
minv(0, 2) = (m(0, 1) * m(1, 2) - m(0, 2) * m(1, 1)) * invdet;
minv(1, 0) = (m(1, 2) * m(2, 0) - m(1, 0) * m(2, 2)) * invdet;
minv(1, 1) = (m(0, 0) * m(2, 2) - m(0, 2) * m(2, 0)) * invdet;
minv(1, 2) = (m(1, 0) * m(0, 2) - m(0, 0) * m(1, 2)) * invdet;
minv(2, 0) = (m(1, 0) * m(2, 1) - m(2, 0) * m(1, 1)) * invdet;
minv(2, 1) = (m(2, 0) * m(0, 1) - m(0, 0) * m(2, 1)) * invdet;
minv(2, 2) = (m(0, 0) * m(1, 1) - m(1, 0) * m(0, 1)) * invdet;
回答by batty
This piece of code computes the transposed inverseof the matrix A:
这段代码计算矩阵 A的转置逆矩阵：
double determinant = +A(0,0)*(A(1,1)*A(2,2)-A(2,1)*A(1,2))
-A(0,1)*(A(1,0)*A(2,2)-A(1,2)*A(2,0))
+A(0,2)*(A(1,0)*A(2,1)-A(1,1)*A(2,0));
double invdet = 1/determinant;
result(0,0) = (A(1,1)*A(2,2)-A(2,1)*A(1,2))*invdet;
result(1,0) = -(A(0,1)*A(2,2)-A(0,2)*A(2,1))*invdet;
result(2,0) = (A(0,1)*A(1,2)-A(0,2)*A(1,1))*invdet;
result(0,1) = -(A(1,0)*A(2,2)-A(1,2)*A(2,0))*invdet;
result(1,1) = (A(0,0)*A(2,2)-A(0,2)*A(2,0))*invdet;
result(2,1) = -(A(0,0)*A(1,2)-A(1,0)*A(0,2))*invdet;
result(0,2) = (A(1,0)*A(2,1)-A(2,0)*A(1,1))*invdet;
result(1,2) = -(A(0,0)*A(2,1)-A(2,0)*A(0,1))*invdet;
result(2,2) = (A(0,0)*A(1,1)-A(1,0)*A(0,1))*invdet;
Though the question stipulated non-singular matrices, you might still want to check if determinant equals zero (or very near zero) and flag it in some way to be safe.
尽管该问题规定了非奇异矩阵，但您可能仍想检查行列式是否等于零（或非常接近于零）并以某种方式对其进行标记以确保安全。
回答by Mr.Ree
With all due respect to our unknown (yahoo) poster, I look at code like that and just die a little inside. Alphabet soup is just so insanely difficult to debug. A single typo anywhere in there can really ruin your whole day. Sadly, this particular example lacked variables with underscores. It's so much more fun when we have a_b-c_d*e_f-g_h. Especially when using a font where _ and - have the same pixel length.
恕我直言，我们不知名（雅虎）海报，我看着这样的代码，只是在里面死了一点。Alphabet 汤实在是太难调试了。任何地方的一个错字都会毁了你的一整天。遗憾的是，这个特殊的例子缺少带下划线的变量。当我们有 a_b-c_d*e_f-g_h 时会更有趣。特别是在使用 _ 和 - 具有相同像素长度的字体时。
Taking up Suvesh Pratapa on his suggestion, I note:
接受 Suvesh Pratapa 的建议，我注意到：
Given 3x3 matrix:
y0x0 y0x1 y0x2
y1x0 y1x1 y1x2
y2x0 y2x1 y2x2
Declared as double matrix [/*Y=*/3] [/*X=*/3];
(A) When taking a minor of a 3x3 array, we have 4 values of interest. The lower X/Y index is always 0 or 1. The higher X/Y index is always 1 or 2. Always! Therefore:
(A) 当取一个 3x3 数组的次要时，我们有 4 个感兴趣的值。较低的 X/Y 指数始终为 0 或 1。较高的 X/Y 指数始终为 1 或 2。始终！所以：
double determinantOfMinor( int theRowHeightY,
int theColumnWidthX,
const double theMatrix [/*Y=*/3] [/*X=*/3] )
{
int x1 = theColumnWidthX == 0 ? 1 : 0; /* always either 0 or 1 */
int x2 = theColumnWidthX == 2 ? 1 : 2; /* always either 1 or 2 */
int y1 = theRowHeightY == 0 ? 1 : 0; /* always either 0 or 1 */
int y2 = theRowHeightY == 2 ? 1 : 2; /* always either 1 or 2 */
return ( theMatrix [y1] [x1] * theMatrix [y2] [x2] )
- ( theMatrix [y1] [x2] * theMatrix [y2] [x1] );
}
(B) Determinant is now: (Note the minus sign!)
(B) 行列式现在是：（注意减号！）
double determinant( const double theMatrix [/*Y=*/3] [/*X=*/3] )
{
return ( theMatrix [0] [0] * determinantOfMinor( 0, 0, theMatrix ) )
- ( theMatrix [0] [1] * determinantOfMinor( 0, 1, theMatrix ) )
+ ( theMatrix [0] [2] * determinantOfMinor( 0, 2, theMatrix ) );
}
(C) And the inverse is now:
(C) 而现在的倒数是：
bool inverse( const double theMatrix [/*Y=*/3] [/*X=*/3],
double theOutput [/*Y=*/3] [/*X=*/3] )
{
double det = determinant( theMatrix );
/* Arbitrary for now. This should be something nicer... */
if ( ABS(det) < 1e-2 )
{
memset( theOutput, 0, sizeof theOutput );
return false;
}
double oneOverDeterminant = 1.0 / det;
for ( int y = 0; y < 3; y ++ )
for ( int x = 0; x < 3; x ++ )
{
/* Rule is inverse = 1/det * minor of the TRANSPOSE matrix. *
* Note (y,x) becomes (x,y) INTENTIONALLY here! */
theOutput [y] [x]
= determinantOfMinor( x, y, theMatrix ) * oneOverDeterminant;
/* (y0,x1) (y1,x0) (y1,x2) and (y2,x1) all need to be negated. */
if( 1 == ((x + y) % 2) )
theOutput [y] [x] = - theOutput [y] [x];
}
return true;
}
And round it out with a little lower-quality testing code:
并用一些质量较低的测试代码来完善它：
void printMatrix( const double theMatrix [/*Y=*/3] [/*X=*/3] )
{
for ( int y = 0; y < 3; y ++ )
{
cout << "[ ";
for ( int x = 0; x < 3; x ++ )
cout << theMatrix [y] [x] << " ";
cout << "]" << endl;
}
cout << endl;
}
void matrixMultiply( const double theMatrixA [/*Y=*/3] [/*X=*/3],
const double theMatrixB [/*Y=*/3] [/*X=*/3],
double theOutput [/*Y=*/3] [/*X=*/3] )
{
for ( int y = 0; y < 3; y ++ )
for ( int x = 0; x < 3; x ++ )
{
theOutput [y] [x] = 0;
for ( int i = 0; i < 3; i ++ )
theOutput [y] [x] += theMatrixA [y] [i] * theMatrixB [i] [x];
}
}
int
main(int argc, char **argv)
{
if ( argc > 1 )
SRANDOM( atoi( argv[1] ) );
double m[3][3] = { { RANDOM_D(0,1e3), RANDOM_D(0,1e3), RANDOM_D(0,1e3) },
{ RANDOM_D(0,1e3), RANDOM_D(0,1e3), RANDOM_D(0,1e3) },
{ RANDOM_D(0,1e3), RANDOM_D(0,1e3), RANDOM_D(0,1e3) } };
double o[3][3], mm[3][3];
if ( argc <= 2 )
cout << fixed << setprecision(3);
printMatrix(m);
cout << endl << endl;
SHOW( determinant(m) );
cout << endl << endl;
BOUT( inverse(m, o) );
printMatrix(m);
printMatrix(o);
cout << endl << endl;
matrixMultiply (m, o, mm );
printMatrix(m);
printMatrix(o);
printMatrix(mm);
cout << endl << endl;
}
Afterthought:
事后思考：
You may also want to detect very large determinants as round-off errors will affect your accuracy!
您可能还想检测非常大的决定因素，因为舍入误差会影响您的准确性！
回答by epatel
A rather nice (I think) header file containing macros for most 2x2, 3x3 and 4x4 matrix operations has been available with most OpenGL toolkits. Not as standard but I've seen it at various places.
大多数 OpenGL 工具包都提供了一个相当不错的（我认为）头文件，其中包含大多数 2x2、3x3 和 4x4 矩阵运算的宏。不是标准的，但我在很多地方都看到过。
You can check it out here. At the end of it you will find both inverse of 2x2, 3x3 and 4x4.
你可以在这里查看。最后你会发现 2x2、3x3 和 4x4 的倒数。
回答by Michael Anderson
Don't try to do this yourself if you're serious about getting edge cases right. So while they many naive/simple methods are theoretically exact, they can have nasty numerical behavior for nearly singular matrices. In particular you can get cancelation/round-off errors that cause you to get arbitrarily bad results.
如果您真的想正确处理边缘情况，请不要尝试自己执行此操作。因此，尽管许多幼稚/简单的方法在理论上是准确的，但对于几乎奇异的矩阵，它们可能具有令人讨厌的数值行为。特别是您可能会遇到取消/舍入错误，从而导致您获得任意糟糕的结果。
A "correct" way is Gaussian elimination with row and column pivoting so that you're always dividing by the largest remaining numerical value. (This is also stable for NxN matrices.). Note that row pivoting alone doesn't catch all the bad cases.
“正确”的方法是使用行和列旋转的高斯消元法，以便您始终除以最大的剩余数值。（这对于 NxN 矩阵也是稳定的。）。请注意，单独行旋转并不能捕获所有坏情况。
However IMO implementing this right and fast is not worth your time - use a well tested library and there are a heap of header only ones.
然而，IMO 正确且快速地实现这一点并不值得您花时间 - 使用经过良好测试的库，并且只有一堆头文件。
回答by moldovean
I have just created a QMatrix class. It uses the built in vector > container. QMatrix.hIt uses the Jordan-Gauss method to compute the inverse of a square matrix.
我刚刚创建了一个 QMatrix 类。它使用内置的 vector > 容器。QMatrix.h它使用 Jordan-Gauss 方法来计算方阵的逆矩阵。
You can use it as follows:
您可以按如下方式使用它：
#include "QMatrix.h"
#include <iostream>
int main(){
QMatrix<double> A(3,3,true);
QMatrix<double> Result = A.inverse()*A; //should give the idendity matrix
std::cout<<A.inverse()<<std::endl;
std::cout<<Result<<std::endl; // for checking
return 0;
}
The inverse function is implemented as follows:
反函数实现如下：
Given a class with the following fields:
给定一个具有以下字段的类：
template<class T> class QMatrix{
public:
int rows, cols;
std::vector<std::vector<T> > A;
the inverse() function:
inverse() 函数：
template<class T>
QMatrix<T> QMatrix<T>:: inverse(){
Identity<T> Id(rows); //the Identity Matrix as a subclass of QMatrix.
QMatrix<T> Result = *this; // making a copy and transforming it to the Identity matrix
T epsilon = 0.000001;
for(int i=0;i<rows;++i){
//check if Result(i,i)==0, if true, switch the row with another
for(int j=i;j<rows;++j){
if(std::abs(Result(j,j))<epsilon) { //uses Overloading()(int int) to extract element from Result Matrix
Result.replace_rows(i,j+1); //switches rows i with j+1
}
else break;
}
// main part, making a triangular matrix
Id(i)=Id(i)*(1.0/Result(i,i));
Result(i)=Result(i)*(1.0/Result(i,i)); // Using overloading ()(int) to get a row form the matrix
for(int j=i+1;j<rows;++j){
T temp = Result(j,i);
Result(j) = Result(j) - Result(i)*temp;
Id(j) = Id(j) - Id(i)*temp; //doing the same operations to the identity matrix
Result(j,i)=0; //not necessary, but looks nicer than 10^-15
}
}
// solving a triangular matrix
for(int i=rows-1;i>0;--i){
for(int j=i-1;j>=0;--j){
T temp = Result(j,i);
Id(j) = Id(j) - temp*Id(i);
Result(j)=Result(j)-temp*Result(i);
}
}
return Id;
}
回答by Larry Gritz
回答by aqeel
# include <conio.h>
# include<iostream.h>
const int size = 9;
int main()
{
char ch;
do
{
clrscr();
int i, j, x, y, z, det, a[size], b[size];
cout << " **** MATRIX OF 3x3 ORDER ****"
<< endl
<< endl
<< endl;
for (i = 0; i <= size; i++)
a[i]=0;
for (i = 0; i < size; i++)
{
cout << "Enter "
<< i + 1
<< " element of matrix=";
cin >> a[i];
cout << endl
<<endl;
}
clrscr();
cout << "your entered matrix is "
<< endl
<<endl;
for (i = 0; i < size; i += 3)
cout << a[i]
<< " "
<< a[i+1]
<< " "
<< a[i+2]
<< endl
<<endl;
cout << "Transpose of given matrix is"
<< endl
<< endl;
for (i = 0; i < 3; i++)
cout << a[i]
<< " "
<< a[i+3]
<< " "
<< a[i+6]
<< endl
<< endl;
cout << "Determinent of given matrix is = ";
x = a[0] * (a[4] * a[8] -a [5] * a[7]);
y = a[1] * (a[3] * a[8] -a [5] * a[6]);
z = a[2] * (a[3] * a[7] -a [4] * a[6]);
det = x - y + z;
cout << det
<< endl
<< endl
<< endl
<< endl;
if (det == 0)
{
cout << "As Determinent=0 so it is singular matrix and its inverse cannot exist"
<< endl
<< endl;
goto quit;
}
b[0] = a[4] * a[8] - a[5] * a[7];
b[1] = a[5] * a[6] - a[3] * a[8];
b[2] = a[3] * a[7] - a[4] * a[6];
b[3] = a[2] * a[7] - a[1] * a[8];
b[4] = a[0] * a[8] - a[2] * a[6];
b[5] = a[1] * a[6] - a[0] * a[7];
b[6] = a[1] * a[5] - a[2] * a[4];
b[7] = a[2] * a[3] - a[0] * a[5];
b[8] = a[0] * a[4] - a[1] * a[3];
cout << "Adjoint of given matrix is"
<< endl
<< endl;
for (i = 0; i < 3; i++)
{
cout << b[i]
<< " "
<< b[i+3]
<< " "
<< b[i+6]
<< endl
<<endl;
}
cout << endl
<<endl;
cout << "Inverse of given matrix is "
<< endl
<< endl
<< endl;
for (i = 0; i < 3; i++)
{
cout << b[i]
<< "/"
<< det
<< " "
<< b[i+3]
<< "/"
<< det
<< " "
<< b[i+6]
<< "/"
<< det
<< endl
<<endl;
}
quit:
cout << endl
<< endl;
cout << "Do You want to continue this again press (y/yes,n/no)";
cin >> ch;
cout << endl
<< endl;
} /* end do */
while (ch == 'y');
getch ();
return 0;
}
回答by Matthew
#include <iostream>
using namespace std;
int main()
{
double A11, A12, A13;
double A21, A22, A23;
double A31, A32, A33;
double B11, B12, B13;
double B21, B22, B23;
double B31, B32, B33;
cout << "Enter all number from left to right, from top to bottom, and press enter after every number: ";
cin >> A11;
cin >> A12;
cin >> A13;
cin >> A21;
cin >> A22;
cin >> A23;
cin >> A31;
cin >> A32;
cin >> A33;
B11 = 1 / ((A22 * A33) - (A23 * A32));
B12 = 1 / ((A13 * A32) - (A12 * A33));
B13 = 1 / ((A12 * A23) - (A13 * A22));
B21 = 1 / ((A23 * A31) - (A21 * A33));
B22 = 1 / ((A11 * A33) - (A13 * A31));
B23 = 1 / ((A13 * A21) - (A11 * A23));
B31 = 1 / ((A21 * A32) - (A22 * A31));
B32 = 1 / ((A12 * A31) - (A11 * A32));
B33 = 1 / ((A11 * A22) - (A12 * A21));
cout << B11 << "\t" << B12 << "\t" << B13 << endl;
cout << B21 << "\t" << B22 << "\t" << B23 << endl;
cout << B31 << "\t" << B32 << "\t" << B33 << endl;
return 0;
}
