Я нашел некоторый код на С++ для определения определителя матрицы, от 4x4 до 8x8. Он работает нормально, но для моего проекта нужны матрицы размером 18x18 или более, а код слишком медленный. Код рекурсивный, но рекурсия правильная концепция для работы с матрицей 18x18? Как еще я могу найти определитель?
Как найти определитель большой матрицы
Ответ 1
Я предполагаю, что вы используете наивный метод расширения формулы Лапласа. Если вы хотите получить скорость, вы можете разложить матрицу M
с помощью LU-декомпозиции (в две матрицы с нижней и верхней диагональю) которую вы можете достичь с помощью модифицированного исключения Гаусса-Джордана в 2*n^3/3 FLOPS
, а затем вычислить определитель как:
det(M) = det(L) * det(U)
, что для треугольных матриц является просто произведением записей в их диагонали.
Этот процесс будет выполняться быстрее, чем O(n!)
.
Изменить: вы также можете использовать метод Crout, который широко внедряется.
Ответ 2
Ну, не многие из нас, работающие в этой области, рассматривали бы 18x18 как большую матрицу, и почти любой метод, который вы выберете, должен быть достаточно быстрым на любом современном компьютере. Также многие из нас не решат матричные вопросы с помощью рекурсивных алгоритмов, гораздо более вероятно, будут использовать итеративные, но это может быть отражением того факта, что многие люди, работающие над матричными проблемами, являются учеными и инженерами, а не учеными-компьютерщиками.
Я предлагаю вам посмотреть на Numerical Recipes в С++. Не обязательно лучший код, который вы найдете, но это текст для изучения и обучения. Для лучших кодов BOOST имеет хорошую репутацию, и всегда есть BLAS и такие вещи, как библиотека ядра Intel Maths или библиотека ядра AMD Core Maths. Я думаю, что все они имеют реализации процедур определения детерминанта, которые очень быстро решат матрицу 18x18.
Ответ 3
Поскольку я не могу комментировать, я хочу добавить это: разложение Холески (или его вариант, LDL T L - нижняя треугольная матрица и D - диагональная матрица) можно использовать для проверьте, что симметричная матрица положительно/отрицательно определена: если она положительно определена, элементы D все положительны, а разложение Холецкого завершится успешно, не беря квадратный корень отрицательного числа. Если матрица отрицательно определена, элементы D все отрицательны, сама матрица не будет иметь декомпозицию Холецкого, но ее отрицательный результат.
"Вычисление определителя треугольной матрицы прост: умножьте диагональные элементы, так как кофакторы недиагональных членов равны 0. Используя разложение LU, это еще упрощает это, так как L - единичная нижнетреугольная матрица, т.е. его диагональные элементы - все 1, в большинстве реализаций. Поэтому вам часто приходится вычислять детерминант U.
- Вы забыли здесь принять во внимание, что все практические реализации гауссовой элиминации используют (частичную) поворот для дополнительной численной устойчивости; поэтому ваше описание является неполным; один подсчитывает количество сводных свопов, выполненных во время фазы декомпозиции, и после умножения всех диагональных элементов U, этот продукт следует сбрасывать, если количество свопов нечетное.
Что касается кода, NR не является бесплатным; Я предлагаю взглянуть на LAPACK/CLAPACK/LAPACK ++ @http://www.netlib.org/. Для справки я могу сделать не лучше, чем указать на "Матричные вычисления" Голуб и Ван Лоан.
Ответ 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();
}