Ho trovato un codice C++ per trovare il determinante della matrice, da 4x4 a 8x8. Funziona bene, ma il mio progetto ha bisogno di matrici di 18x18 o più e il codice è troppo lento. Il codice è ricorsivo, ma la ricorsione è il concetto giusto per trattare una matrice 18x18? In quale altro modo posso trovare il determinante?Come trovare il determinante della matrice grande
Suppongo che tu stia utilizzando il metodo naive di espansione Laplace's formula. Se si vuole guadagnare in velocità, è possibile scomporre tua matrice M utilizzando LU decomposition (in due matrici inferiore e superiore diagonale), che si può ottenere con un'eliminazione di Gauss-Jordan modificato in 2*n^3/3 FLOPS e quindi calcolare il determinante come:
det(M) = det(L) * det(U), che per le matrici triangolari è solo il prodotto delle voci nella loro diagonale.
Questo processo sarà ancora più veloce di O(n!).
Modifica: è inoltre possibile utilizzare Crout's method, che è ampiamente implementato.
Se la matrice è simmetrica positiva-definita, la fattorizzazione di Cholesky è il metodo più veloce e numericamente migliore. – Paul
@Paul: È vero, ma non sono sicuro di come verificare che una matrice sia definita positiva, quindi l'ho saltata. –
Calcolare il determinante di una matrice triangolare è semplice: moltiplica gli elementi diagonali, come i cofattori dei termini fuori diagonale sono 0. L'uso di una decomposizione LU lo semplifica ulteriormente, poiché L è un'unità, matrice triangolare inferiore, cioè i suoi elementi diagonali sono tutti 1, nella maggior parte delle implementazioni. Perciò, spesso devi calcolare solo il determinante di U. – rcollyer
Bene, non molti di noi che lavorano nel campo considererebbero 18x18 come una matrice di grandi dimensioni e quasi tutte le tecniche che si sceglie dovrebbero essere abbastanza veloci su qualsiasi computer moderno. Né molti di noi affrontano questioni a matrice con algoritmi ricorsivi, molto più propensi a usare quelli iterativi, ma ciò potrebbe riflettere il fatto che molte persone che lavorano su problemi con la matrice sono scienziati e ingegneri, non scienziati informatici.
Suggerisco di dare un'occhiata alle Ricette numeriche in C++. Non necessariamente il miglior codice che troverai, ma è un testo da studiare e da cui imparare. Per codici migliori, BOOST ha una buona reputazione e c'è sempre BLAS e cose come Intel Maths Kernel Library o AMD Core Maths Library. Penso che tutti questi abbiano implementazioni di routine di determinazione dei determinanti che affronteranno una matrice 18x18 molto rapidamente.
Poiché non posso commentare, desidero aggiungere questo: la decomposizione di Cholesky (o sua variante, LDL T, L un'unità inferiore matrice triangolare e D una matrice diagonale) può essere utilizzato per verificare se un simmetrica la matrice è definita positiva/negativa: se è definita positiva, gli elementi di D sono tutti positivi e la decomposizione di Cholesky terminerà correttamente senza prendere la radice quadrata di un numero negativo. Se la matrice è definita come negativa, gli elementi di D sono tutti negativi, la matrice stessa non avrà una decomposizione di Cholesky, ma il negativo di essa sarebbe.
"Calcolare il determinante di una matrice triangolare è semplice: moltiplica gli elementi diagonali, come i cofattori dei termini fuori diagonale sono 0. L'uso di una decomposizione LU lo semplifica ulteriormente, poiché L è un'unità, matrice triangolare inferiore, cioè i suoi elementi diagonali sono tutti 1, nella maggior parte delle implementazioni, quindi spesso devi solo calcolare il determinante di U. "
Come per il codice, NR non è libero; Suggerisco di dare un'occhiata a LAPACK/CLAPACK/LAPACK ++ @http://www.netlib.org/ invece. Come riferimento, non posso fare di meglio che indirizzarti a "Matrix Computations" di Golub e Van Loan.
Penso che questo potrebbe funzionare.
L'ho scritto quando ho studiato corso di analisi numerica.
questo non è solo determinanti, ma altre funzioni relative alla matrice
In primo luogo, copiare e salvare il codice come '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;
}
Nel mio metodo, ho convertire il matrice alla matrice triangolare superiore utilizzando il funzionamento di riga elementare
E il determinante è il prodotto degli elementi diagonali.
Ecco il codice di esempio
#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();
}
codice Python La funzione det_recursive lavora per una matrice quadrata di qualsiasi dimensione. Tuttavia, dal momento che sta usando il metodo ricorsive naive di espansione della formula di Laplace, è molto lento per le matrici di grandi dimensioni.
Un'altra tecnica consiste nel trasformare la matrice in una forma triangolare superiore utilizzando la tecnica di eliminazione di gauss. Quindi il determinante della matrice è semplicemente il prodotto di elementi diagonali della forma trasformata triangolare della matrice originale.
Fondamentalmente numpy è il metodo più veloce, ma internamente utilizza una sorta di metodo di trasformazione della matrice lineare simile a quello che elimina la gauss. Tuttavia, non sono sicuro di cosa sia esattamente!
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)
@Pegah Per rispondere a una domanda, scorrere fino a "La tua risposta". Ci piace mantenere separate le domande e le risposte. – phihag