大きな行列の行列式を見つける方法

Question

4x4から8x8の行列の行列式を見つけるためのC++コードが見つかりました。それは正常に動作しますが、私のプロジェクトは18x18以上の行列を必要とし、コードが遅すぎます。コードは再帰的ですが、再帰は18x18行列を扱う正しい概念ですか？どうやって行列式を見つけることができますか？

Answer 1

ナイーブメソッドを展開中のLaplace's formulaとしているとします。

：あなたはスピードで獲得したい場合は、あなたが 2*n^3/3 FLOPSで修正ガウス-Jordanの消去を達成し、次にような決定を計算することができ LU decomposition（2小文字と上対角行列への）を使用して行列 Mを分解することができます

det(M) = det(L) * det(U)であり、これは三角行列の場合、対角線上のエントリの積である。

このプロセスはまだO(n!)より速くなります。

編集：Crout's methodも広く実装されています。

Answer 2

現場で働いている私たちの多くは、18x18を大規模なマトリックスとみなしておらず、あなたが選んだほとんどすべての技術は現代のコンピュータでは十分速くなければなりません。私たちの多くは、反復アルゴリズムを使用して行列質問に取り組むことはありませんが、反復アルゴリズムを使用する可能性は非常に高いですが、行列問題を扱う多くの人々がコンピュータ科学者ではなく科学者やエンジニアであるという事実を反映しています。

C++の数値レシピをお勧めします。必ずしも最良のコードではありませんが、勉強や学習のためのテキストです。より良いコードのためには、BOOSTは評判がよく、BLASや、Intel Maths Kernel LibraryやAMD Core Maths Libraryなどがあります。私はこれらのすべてが非常に迅速に18x18行列に取り組む行列式決定ルーチンの実装を持っていると思います。

Answer 3

コレスキ分解（またはその変形、LDL ^T、Lは単位下三角行列、Dは対角行列）を使用して、対称行列は正または負の定理である：正定値ならば、Dの要素はすべて正であり、コレスキー分解は負の数の平方根を取ることなく首尾よく終了する。行列が負定値である場合、Dの要素はすべて負であり、行列自体はコレスキー分解を持たないが、負の値をとる。

"三角行列の行列式を計算するのは簡単です：対角要素の補因子が0であるため、対角要素を掛けます。LU分解を使用すると、Lが単位、下三角行列、つまり、その対角要素はほとんどの実装では1です。そのため、Uの行列式を計算するだけでよい場合があります。

• ガウス消去の実用的な実装では、追加の数値安定性のために（部分）ピボットを利用することを考慮してここを忘れています。あなたの説明は不完全です。分解フェーズ中に実行された行スワップの数をカウントし、Uのすべての対角要素を掛け合わせた後、スワップの数が奇数の場合、この積は否定されるべきである。

NRはフリーではありません。代わりにLAPACK/CLAPACK/LAPACK ++ @http://www.netlib.org/を見てみることをお勧めします。参考までに、私はGolubとVan Loanの "Matrix Computations"を指摘するよりもうまくやっていけません。

Answer 4

これはおそらく動作すると思います。

数値解析コースを勉強したときに書きました。

これは私が変換し、私の方法では「Matrix.h」

//Title: Matrix Header File 
//Writer: Say OL 
//This is a beginner code not an expert one 
//No responsibilty for any errors 
//Use for your own risk 
using namespace std; 
int row,col,Row,Col; 
double Coefficient; 
//Input Matrix 
void Input(double Matrix[9][9],int Row,int Col) 
{ 
    for(row=1;row<=Row;row++) 
     for(col=1;col<=Col;col++) 
     { 
      cout<<"e["<<row<<"]["<<col<<"]="; 
      cin>>Matrix[row][col]; 
     } 
} 
//Output Matrix 
void Output(double Matrix[9][9],int Row,int Col) 
{ 
    for(row=1;row<=Row;row++) 
    { 
     for(col=1;col<=Col;col++) 
      cout<<Matrix[row][col]<<"\t"; 
     cout<<endl; 
    } 
} 
//Copy Pointer to Matrix 
void CopyPointer(double (*Pointer)[9],double Matrix[9][9],int Row,int Col) 
{ 
    for(row=1;row<=Row;row++) 
     for(col=1;col<=Col;col++) 
      Matrix[row][col]=Pointer[row][col]; 
} 
//Copy Matrix to Matrix 
void CopyMatrix(double MatrixInput[9][9],double MatrixTarget[9][9],int Row,int Col) 
{ 
    for(row=1;row<=Row;row++) 
     for(col=1;col<=Col;col++) 
      MatrixTarget[row][col]=MatrixInput[row][col]; 
} 
//Transpose of Matrix 
double MatrixTran[9][9]; 
double (*(Transpose)(double MatrixInput[9][9],int Row,int Col))[9] 
{ 
    for(row=1;row<=Row;row++) 
     for(col=1;col<=Col;col++) 
      MatrixTran[col][row]=MatrixInput[row][col]; 
    return MatrixTran; 
} 
//Matrix Addition 
double MatrixAdd[9][9]; 
double (*(Addition)(double MatrixA[9][9],double MatrixB[9][9],int Row,int Col))[9] 
{ 
    for(row=1;row<=Row;row++) 
     for(col=1;col<=Col;col++) 
      MatrixAdd[row][col]=MatrixA[row][col]+MatrixB[row][col]; 
    return MatrixAdd; 
} 
//Matrix Subtraction 
double MatrixSub[9][9]; 
double (*(Subtraction)(double MatrixA[9][9],double MatrixB[9][9],int Row,int Col))[9] 
{ 
    for(row=1;row<=Row;row++) 
     for(col=1;col<=Col;col++) 
      MatrixSub[row][col]=MatrixA[row][col]-MatrixB[row][col]; 
    return MatrixSub; 
} 
//Matrix Multiplication 
int mRow,nCol,pCol,kcol; 
double MatrixMult[9][9]; 
double (*(Multiplication)(double MatrixA[9][9],double MatrixB[9][9],int mRow,int nCol,int pCol))[9] 
{ 
    for(row=1;row<=mRow;row++) 
     for(col=1;col<=pCol;col++) 
     { 
      MatrixMult[row][col]=0.0; 
      for(kcol=1;kcol<=nCol;kcol++) 
       MatrixMult[row][col]+=MatrixA[row][kcol]*MatrixB[kcol][col]; 
     } 
    return MatrixMult; 
} 
//Interchange Two Rows 
double RowTemp[9][9]; 
double MatrixInter[9][9]; 
double (*(InterchangeRow)(double MatrixInput[9][9],int Row,int Col,int iRow,int jRow))[9] 
{ 
    CopyMatrix(MatrixInput,MatrixInter,Row,Col); 
    for(col=1;col<=Col;col++) 
    { 
     RowTemp[iRow][col]=MatrixInter[iRow][col]; 
     MatrixInter[iRow][col]=MatrixInter[jRow][col]; 
     MatrixInter[jRow][col]=RowTemp[iRow][col]; 
    } 
    return MatrixInter; 
} 
//Pivote Downward 
double MatrixDown[9][9]; 
double (*(PivoteDown)(double MatrixInput[9][9],int Row,int Col,int tRow,int tCol))[9] 
{ 
    CopyMatrix(MatrixInput,MatrixDown,Row,Col); 
    Coefficient=MatrixDown[tRow][tCol]; 
    if(Coefficient!=1.0) 
     for(col=1;col<=Col;col++) 
      MatrixDown[tRow][col]/=Coefficient; 
    if(tRow<Row) 
     for(row=tRow+1;row<=Row;row++) 
     { 
      Coefficient=MatrixDown[row][tCol]; 
      for(col=1;col<=Col;col++) 
       MatrixDown[row][col]-=Coefficient*MatrixDown[tRow][col]; 
     } 
return MatrixDown; 
} 
//Pivote Upward 
double MatrixUp[9][9]; 
double (*(PivoteUp)(double MatrixInput[9][9],int Row,int Col,int tRow,int tCol))[9] 
{ 
    CopyMatrix(MatrixInput,MatrixUp,Row,Col); 
    Coefficient=MatrixUp[tRow][tCol]; 
    if(Coefficient!=1.0) 
     for(col=1;col<=Col;col++) 
      MatrixUp[tRow][col]/=Coefficient; 
    if(tRow>1) 
     for(row=tRow-1;row>=1;row--) 
     { 
      Coefficient=MatrixUp[row][tCol]; 
      for(col=1;col<=Col;col++) 
       MatrixUp[row][col]-=Coefficient*MatrixUp[tRow][col]; 
     } 
    return MatrixUp; 
} 
//Pivote in Determinant 
double MatrixPiv[9][9]; 
double (*(Pivote)(double MatrixInput[9][9],int Dim,int pTarget))[9] 
{ 
    CopyMatrix(MatrixInput,MatrixPiv,Dim,Dim); 
    for(row=pTarget+1;row<=Dim;row++) 
    { 
     Coefficient=MatrixPiv[row][pTarget]/MatrixPiv[pTarget][pTarget]; 
     for(col=1;col<=Dim;col++) 
     { 
      MatrixPiv[row][col]-=Coefficient*MatrixPiv[pTarget][col]; 
     } 
    } 
    return MatrixPiv; 
} 
//Determinant of Square Matrix 
int dCounter,dRow; 
double Det; 
double MatrixDet[9][9]; 
double Determinant(double MatrixInput[9][9],int Dim) 
{ 
    CopyMatrix(MatrixInput,MatrixDet,Dim,Dim); 
    Det=1.0; 
    if(Dim>1) 
    { 
     for(dRow=1;dRow<Dim;dRow++) 
     { 
      dCounter=dRow; 
      while((MatrixDet[dRow][dRow]==0.0)&(dCounter<=Dim)) 
      { 
       dCounter++; 
       Det*=-1.0; 
       CopyPointer(InterchangeRow(MatrixDet,Dim,Dim,dRow,dCounter),MatrixDet,Dim,Dim); 
      } 
      if(MatrixDet[dRow][dRow]==0) 
      { 
       Det=0.0; 
       break; 
      } 
      else 
      { 
       Det*=MatrixDet[dRow][dRow]; 
       CopyPointer(Pivote(MatrixDet,Dim,dRow),MatrixDet,Dim,Dim); 
      } 
     } 
     Det*=MatrixDet[Dim][Dim]; 
    } 
    else Det=MatrixDet[1][1]; 
    return Det; 
} 
//Matrix Identity 
double MatrixIdent[9][9]; 
double (*(Identity)(int Dim))[9] 
{ 
    for(row=1;row<=Dim;row++) 
     for(col=1;col<=Dim;col++) 
      if(row==col) 
       MatrixIdent[row][col]=1.0; 
      else 
       MatrixIdent[row][col]=0.0; 
    return MatrixIdent; 
} 
//Join Matrix to be Augmented Matrix 
double MatrixJoin[9][9]; 
double (*(JoinMatrix)(double MatrixA[9][9],double MatrixB[9][9],int Row,int ColA,int ColB))[9] 
{ 
    Col=ColA+ColB; 
    for(row=1;row<=Row;row++) 
     for(col=1;col<=Col;col++) 
      if(col<=ColA) 
       MatrixJoin[row][col]=MatrixA[row][col]; 
      else 
       MatrixJoin[row][col]=MatrixB[row][col-ColA]; 
    return MatrixJoin; 
} 
//Inverse of Matrix 
double (*Pointer)[9]; 
double IdentMatrix[9][9]; 
int Counter; 
double MatrixAug[9][9]; 
double MatrixInv[9][9]; 
double (*(Inverse)(double MatrixInput[9][9],int Dim))[9] 
{ 
    Row=Dim; 
    Col=Dim+Dim; 
    Pointer=Identity(Dim); 
    CopyPointer(Pointer,IdentMatrix,Dim,Dim); 
    Pointer=JoinMatrix(MatrixInput,IdentMatrix,Dim,Dim,Dim); 
    CopyPointer(Pointer,MatrixAug,Row,Col); 
    for(Counter=1;Counter<=Dim;Counter++) 
    { 
     Pointer=PivoteDown(MatrixAug,Row,Col,Counter,Counter); 
     CopyPointer(Pointer,MatrixAug,Row,Col); 
    } 
    for(Counter=Dim;Counter>1;Counter--) 
    { 
     Pointer=PivoteUp(MatrixAug,Row,Col,Counter,Counter); 
     CopyPointer(Pointer,MatrixAug,Row,Col); 
    } 
    for(row=1;row<=Dim;row++) 
     for(col=1;col<=Dim;col++) 
      MatrixInv[row][col]=MatrixAug[row][col+Dim]; 
    return MatrixInv; 
} 
//Gauss-Jordan Elemination 
double MatrixGJ[9][9]; 
double VectorGJ[9][9]; 
double (*(GaussJordan)(double MatrixInput[9][9],double VectorInput[9][9],int Dim))[9] 
{ 
    Row=Dim; 
    Col=Dim+1; 
    Pointer=JoinMatrix(MatrixInput,VectorInput,Dim,Dim,1); 
    CopyPointer(Pointer,MatrixGJ,Row,Col); 
    for(Counter=1;Counter<=Dim;Counter++) 
    { 
     Pointer=PivoteDown(MatrixGJ,Row,Col,Counter,Counter); 
     CopyPointer(Pointer,MatrixGJ,Row,Col); 
    } 
    for(Counter=Dim;Counter>1;Counter--) 
    { 
     Pointer=PivoteUp(MatrixGJ,Row,Col,Counter,Counter); 
     CopyPointer(Pointer,MatrixGJ,Row,Col); 
    } 
    for(row=1;row<=Dim;row++) 
     for(col=1;col<=1;col++) 
      VectorGJ[row][col]=MatrixGJ[row][col+Dim]; 
    return VectorGJ; 
} 
//Generalized Gauss-Jordan Elemination 
double MatrixGGJ[9][9]; 
double VectorGGJ[9][9]; 
double (*(GeneralizedGaussJordan)(double MatrixInput[9][9],double VectorInput[9][9],int Dim,int vCol))[9] 
{ 
    Row=Dim; 
    Col=Dim+vCol; 
    Pointer=JoinMatrix(MatrixInput,VectorInput,Dim,Dim,vCol); 
    CopyPointer(Pointer,MatrixGGJ,Row,Col); 
    for(Counter=1;Counter<=Dim;Counter++) 
    { 
     Pointer=PivoteDown(MatrixGGJ,Row,Col,Counter,Counter); 
     CopyPointer(Pointer,MatrixGGJ,Row,Col); 
    } 
    for(Counter=Dim;Counter>1;Counter--) 
    { 
     Pointer=PivoteUp(MatrixGGJ,Row,Col,Counter,Counter); 
     CopyPointer(Pointer,MatrixGGJ,Row,Col); 
    } 
    for(row=1;row<=Row;row++) 
     for(col=1;col<=vCol;col++) 
      VectorGGJ[row][col]=MatrixGGJ[row][col+Dim]; 
    return VectorGGJ; 
} 
//Matrix Sparse, Three Diagonal Non-Zero Elements 
double MatrixSpa[9][9]; 
double (*(Sparse)(int Dimension,double FirstElement,double SecondElement,double ThirdElement))[9] 
{ 
    MatrixSpa[1][1]=SecondElement; 
    MatrixSpa[1][2]=ThirdElement; 
    MatrixSpa[Dimension][Dimension-1]=FirstElement; 
    MatrixSpa[Dimension][Dimension]=SecondElement; 
    for(int Counter=2;Counter<Dimension;Counter++) 
    { 
     MatrixSpa[Counter][Counter-1]=FirstElement; 
     MatrixSpa[Counter][Counter]=SecondElement; 
     MatrixSpa[Counter][Counter+1]=ThirdElement; 
    } 
    return MatrixSpa; 
} 

ように、コードをコピーして保存し、

まず行列に関連する決定が、他の機能だけではありません基本行操作を使用して行列から上三角行列に変換する。

そして行列式は対角要素の積である。ここ

サンプルコード

#include<iostream> 
#include<conio.h> 
#include"Matrix.h" 
int Dim; 
double Matrix[9][9]; 
int main() 
{ 
    cout<<"Enter matrix dimension: "; 
    cin>>Dim; 
    cout<<"Enter matrix elements:"<<endl; 
    Input(Matrix,Dim,Dim); 
    cout<<"Your matrix:"<<endl; 
    Output(Matrix,Dim,Dim); 
    cout<<"The determinant: "<<Determinant(Matrix,Dim)<<endl; 
    getch(); 
} 

Answer 5

Pythonコード 任意のサイズの正方行列で働いdet_recursive関数です。しかし、ラプラスの式を拡張する再帰的な単純な方法を使用しているので、大きなサイズの行列の場合は非常に遅いです。

もう1つの手法は、ガウス消去技術を使用して行列を上三角形に変換することです。次に、行列の行列式は、元の行列の三角形に変換された形の対角要素の積である。

基本的にnumpyは最も速いですが、内部的にはガウス消去と同様の線形行列変換方法を使用しています。しかし、私はそれが何であるかは分かりません！

In[1] 
import numpy as np 

In[2] 
mat = np.random.rand(9,9) 
print("numpy asnwer = ", np.linalg.det(mat)) 

Out[2] 
numpy asnwer = 0.016770106020608373 

In[3] 
def det_recursive(A): 
    if A.shape[0] != A.shape[1]: 
     raise ValueError('matrix {} is not Square'.format(A)) 

    sol = 0 
    if A.shape != (1,1): 
     for i in range(A.shape[0]): 
      sol = sol + (-1)**i * A[i, 0] * det_recursive(np.delete(np.delete(A, 0, axis= 1), i, axis= 0)) 
     return sol 
    else: 
     return A[0,0] 
​ 

In[4] 
print("recursive asnwer = ", det_recursive(mat)) 

Out[4] 
recursive asnwer = 0.016770106020608397 

In[5] 
def det_gauss_elimination(a,tol=1.0e-9): 
    """ 
    calculate determinant using gauss-elimination method 
    """ 
    a = np.copy(a) 

    assert(a.shape[0] == a.shape[1]) 
    n = a.shape[0] 

    # Set up scale factors 
    s = np.zeros(n) 

    mult = 0 
    for i in range(n): 
     s[i] = max(np.abs(a[i,:])) # find the max of each row 
    for k in range(0, n-1): #pivot row 
     # Row interchange, if needed 
     p = np.argmax(np.abs(a[k:n,k])/s[k:n]) + k 
     if abs(a[p,k]) < tol: 
      print("Matrix is singular") 
      return 0 
     if p != k: 
      a[[k,p],:] = a[[p, k],:] 
      s[k],s[p] = s[p],s[k] 
      mult = mult + 1 
​ 
     # convert a to upper triangular matrix 
     for i in range(k+1,n): 
      if a[i,k] != 0.0: # skip if a(i,k) is already zero 
       lam = a [i,k]/a[k,k] 
       a[i,k:n] = a[i,k:n] - lam*a[k,k:n] 
​ 
    deter = np.prod(np.diag(a))* (-1)**mult 
    return deter 

In[6] 
print("gauss elimination asnwer = ", det_gauss_elimination(mat)) 

Out[6] 
gauss elimination asnwer = 0.016770106020608383 

In[7] 
print("numpy time") 
%timeit -n3 -r3 np.linalg.det(mat) 
print("\nrecursion time") 
%timeit -n3 -r3 det_recursive(mat) 
print("\ngauss_elimination time") 
%timeit -n3 -r3 det_gauss_elimination(mat) 

Out[7] 
numpy time 
40.8 µs ± 17.2 µs per loop (mean ± std. dev. of 3 runs, 3 loops each) 

recursion time 
10.1 s ± 128 ms per loop (mean ± std. dev. of 3 runs, 3 loops each) 

gauss_elimination time 
472 µs ± 106 µs per loop (mean ± std. dev. of 3 runs, 3 loops each)