3x3 逆行列を計算する最も簡単な方法は何ですか?
おそらくCramerのルールを使用して、非特異行列のトリックを実行する短いコードスニペットを探しています。高度に最適化する必要はありません。スピードよりもシンプルさを好みます。追加のライブラリにリンクしたくありません。
118714 次
13 に答える
55
これは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;
于 2013-08-29T07:18:59.033 に答える
31
このコードは、行列 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;
質問では非特異行列が規定されていますが、行列式がゼロ (または非常にゼロに近い) に等しいかどうかを確認し、安全のために何らかの方法でフラグを立てたい場合があります。
于 2009-06-11T23:17:36.327 に答える
27
自分でコーディングしてみませんか？挑戦してください。:)
3×3 行列の場合
(ソース: wolfram.com )
逆行列は
(ソース: wolfram.com )
行列式 |A| が何であるかを知っていると仮定しています。は。
画像 (c) Wolfram|Alphaと mathworld.wolfram (06-11-09, 22.06)
于 2009-06-11T22:16:49.157 に答える
9
私たちの未知の (yahoo) ポスターに敬意を表して、私はそのようなコードを見て、ただ少し心の中で死んでいます。アルファベットのスープは、デバッグが非常に困難です。どこかで 1 つのタイプミスをすると、1 日が本当に台無しになってしまいます。悲しいことに、この特定の例にはアンダースコア付きの変数がありませんでした。a_b-c_d*e_f-g_h があると、とても楽しいです。特に、_ と - のピクセル長が同じフォントを使用する場合。
Suvesh Pratapa の提案を取り上げて、私は次のように述べています。
Given 3x3 matrix:
y0x0 y0x1 y0x2
y1x0 y1x1 y1x2
y2x0 y2x1 y2x2
Declared as double matrix [/*Y=*/3] [/*X=*/3];
(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) 行列式は次のようになりました: (マイナス記号に注意してください!)
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) 逆は次のようになります。
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;
}
そして、少し品質の低いテスト コードで締めくくります。
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;
}
結果論：
丸め誤差が精度に影響するため、非常に大きな行列式を検出することもできます。
于 2009-06-12T11:04:12.167 に答える
3
エッジケースを正しく取得することに真剣に取り組んでいる場合は、これを自分でやろうとしないでください。したがって、多くの素朴で単純な方法は理論的には正確ですが、ほぼ特異な行列に対して厄介な数値動作をする可能性があります。特に、キャンセル/丸めエラーが発生して、勝手に悪い結果が得られる可能性があります。
「正しい」方法は、行と列のピボットを使用したガウス消去法であり、常に残りの最大の数値で除算します。(これは、NxN 行列に対しても安定しています。) 行のピボットだけでは、すべての悪いケースを捉えることができないことに注意してください。
ただし、IMO がこれを適切かつ迅速に実装するのは時間の価値がありません。十分にテストされたライブラリを使用すると、ヘッダーのみのライブラリがたくさんあります。
于 2013-08-08T06:31:21.447 に答える
3
大部分の 2x2、3x3、および 4x4 行列演算用のマクロを含むかなり優れた (私が思うに) ヘッダー ファイルが、ほとんどの OpenGL ツールキットで利用可能です。定番ではありませんが、いろいろなところで見かけました。
ここで確認できます。最後に、2x2、3x3、および 4x4 の逆数が表示されます。
于 2009-12-30T19:35:17.323 に答える
1
# 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;
}
于 2010-05-11T08:27:31.410 に答える
1
OpenEXRの一部である Ilmbase もお勧めします。これは、テンプレート化された 2,3,4 ベクトルおよびマトリックス ルーチンの優れたセットです。
于 2009-12-30T23:37:51.393 に答える
0
//Function for inverse of the input square matrix 'J' of dimension 'dim':
vector<vector<double > > inverseVec33(vector<vector<double > > J, int dim)
{
//Matrix of Minors
vector<vector<double > > invJ(dim,vector<double > (dim));
for(int i=0; i<dim; i++)
{
for(int j=0; j<dim; j++)
{
invJ[i][j] = (J[(i+1)%dim][(j+1)%dim]*J[(i+2)%dim][(j+2)%dim] -
J[(i+2)%dim][(j+1)%dim]*J[(i+1)%dim][(j+2)%dim]);
}
}
//determinant of the matrix:
double detJ = 0.0;
for(int j=0; j<dim; j++)
{ detJ += J[0][j]*invJ[0][j];}
//Inverse of the given matrix.
vector<vector<double > > invJT(dim,vector<double > (dim));
for(int i=0; i<dim; i++)
{
for(int j=0; j<dim; j++)
{
invJT[i][j] = invJ[j][i]/detJ;
}
}
return invJT;
}
void main()
{
//given matrix:
vector<vector<double > > Jac(3,vector<double > (3));
Jac[0][0] = 1; Jac[0][1] = 2; Jac[0][2] = 6;
Jac[1][0] = -3; Jac[1][1] = 4; Jac[1][2] = 3;
Jac[2][0] = 5; Jac[2][1] = 1; Jac[2][2] = -4;`
//Inverse of the matrix Jac:
vector<vector<double > > JacI(3,vector<double > (3));
//call function and store inverse of J as JacI:
JacI = inverseVec33(Jac,3);
}
于 2014-03-04T13:37:09.020 に答える
0
#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;
}
于 2013-02-19T02:39:25.717 に答える
0
このような問題については、C ++よりもはるかに読みやすいと思うので、先に進んでPythonで書きました。機能順はこちらの動画で手で解決する操作順です。これをインポートして、マトリックスで「print_invert」を呼び出すだけです。
def print_invert (matrix):
i_matrix = invert_matrix (matrix)
for line in i_matrix:
print (line)
return
def invert_matrix (matrix):
determinant = str (determinant_of_3x3 (matrix))
cofactor = make_cofactor (matrix)
trans_matrix = transpose_cofactor (cofactor)
trans_matrix[:] = [[str (element) +'/'+ determinant for element in row] for row in trans_matrix]
return trans_matrix
def determinant_of_3x3 (matrix):
multiplication = 1
neg_multiplication = 1
total = 0
for start_column in range (3):
for row in range (3):
multiplication *= matrix[row][(start_column+row)%3]
neg_multiplication *= matrix[row][(start_column-row)%3]
total += multiplication - neg_multiplication
multiplication = neg_multiplication = 1
if total == 0:
total = 1
return total
def make_cofactor (matrix):
cofactor = [[0,0,0],[0,0,0],[0,0,0]]
matrix_2x2 = [[0,0],[0,0]]
# For each element in matrix...
for row in range (3):
for column in range (3):
# ...make the 2x2 matrix in this inner loop
matrix_2x2 = make_2x2_from_spot_in_3x3 (row, column, matrix)
cofactor[row][column] = determinant_of_2x2 (matrix_2x2)
return flip_signs (cofactor)
def make_2x2_from_spot_in_3x3 (row, column, matrix):
c_count = 0
r_count = 0
matrix_2x2 = [[0,0],[0,0]]
# ...make the 2x2 matrix in this inner loop
for inner_row in range (3):
for inner_column in range (3):
if row is not inner_row and inner_column is not column:
matrix_2x2[r_count % 2][c_count % 2] = matrix[inner_row][inner_column]
c_count += 1
if row is not inner_row:
r_count += 1
return matrix_2x2
def determinant_of_2x2 (matrix):
total = matrix[0][0] * matrix [1][1]
return total - (matrix [1][0] * matrix [0][1])
def flip_signs (cofactor):
sign_pos = True
# For each element in matrix...
for row in range (3):
for column in range (3):
if sign_pos:
sign_pos = False
else:
cofactor[row][column] *= -1
sign_pos = True
return cofactor
def transpose_cofactor (cofactor):
new_cofactor = [[0,0,0],[0,0,0],[0,0,0]]
for row in range (3):
for column in range (3):
new_cofactor[column][row] = cofactor[row][column]
return new_cofactor
于 2013-03-02T12:51:04.353 に答える