如何在c中实现拉普拉斯展开算法



我很难找到让这个算法工作的方法,因为我不知道如何解决问题的中间部分。到目前为止,这是我的代码:

int det(int matrixSize, int matrix[][matrixSize]){
int determinant = 0, matrixValues[matrixSize * matrixSize], matrixFirstRowValues[matrixSize * matrixSize];
for(int i = matrixSize; i > 2; i--){
for(int row = 0; row < matrixSize; row++){
for(int col = 0; col < matrixSize; col++){
matrixFirstRowValues[row + (matrixSize - i)] = matrix[1][col + (matrixSize - i)]; 
//copies the first row values for an array until we have a 2x2 matrix
}
} 
}
//multiply the values in matrix Values by their respective matrix without 
//the row and column of these values times -1^row+col
determinant += (matrix[matrixSize-1][matrixSize-1] * matrix[matrixSize-2][matrixSize-2])
- (matrix[matrixSize-1][matrixSize-2] * matrix[matrixSize-2][matrixSize-1]);
return determinant;
}

作为矩阵,一个大小为matrixSize的二维数组,我对它进行迭代,直到剩下一个2x2矩阵,将第一行的每个值复制到一个新数组中。

这些值必须乘以当我删除该值所在的行和列时留下的矩阵。

这就是拉普拉斯展开的原理。给我带来麻烦的部分是获取那些通过删除行和列而留下的矩阵,因为我希望这对nxn矩阵有效。

然后,在最后,求和到2x2矩阵的det。

如何使用当前设置完成中间部分(注释所在的位置)?

当我删除该值所在的行和列。

您必须与辅因子矩阵相乘,辅因子矩阵的条目是移除第i行和第j列时剩余矩阵的行列式

当然,这是递归算法的设置,因为较大矩阵的行列式是用较小矩阵的行列式表示的:如果A = (a_{ij})是矩阵,那么det(A) = sum j = 1..n: a_{ij} * det(M_{ij}),其中M_{ij}是在去除i固定的i行和j列时从a产生的次矩阵。基本情况是2乘2矩阵,也可能是3乘3矩阵。

出现的问题是,n-by-n矩阵产生大小为(n-1)-by-(n-1)n矩阵M_{ij},每个矩阵产生大小少一的n-1矩阵,依此类推,直到达到基本情况,此时您必须跟踪n!/2矩阵。(很明显,拉普拉斯展开是一种成本相当高的算法,任何基于高斯消去的算法都会更高效。但这只是一个题外话,因为我们正在讨论拉普拉斯展开。)如果以迭代的方式进行,则必须手动进行,递归算法将通过堆栈框架进行隐式记账。

您的方法

让我们来看看您提供的这段代码。我搞不懂你究竟想达到什么目的。以最内层循环中的语句为例,该语句在col:上迭代

matrixFirstRowValues[row + (matrixSize - i)] = matrix[1][col + (matrixSize - i)];

只要col在最里面的循环中发生变化,rowi都是固定的,所以您将从matrix中的第二行(显然)分配和重新分配到matrixFirstRowValues中的同一个条目。不仅如此,您还可以从索引范围(matrixSize-i) .. (2*matrixSize - (i+1))进行分配,该范围超出了列的范围,除非是i == matrixSize,这只是i的第一个值的情况。

正如我之前提到的,最终你不会得到一个2乘2的矩阵,而是n!/2

除第i行和第j列外的复制

查看删除了i行和j列的矩阵,您会得到四个子矩阵(其中一些可能是空的)。将自己限制为沿第一行展开,然后只处理两个子矩阵(其中一些可能是空的)。您可以使用两个循环,一个用于第j列左侧的矩阵,另一个用于右侧的矩阵,或者,如前一个答案中所建议的,选择使用continue跳过第j列来循环循环循环,而不更新目标列索引。如果col标记要删除的当前列(当前行始终为0),则在所有行上迭代r,在所有列上迭代c,并在c == col处将列循环分解为两部分。假设目标矩阵被称为minor,那么它看起来像这样:

// expand along first row
for(col = 0; col < matrix_size; col++) {
// copy into minor matrix, left side then right side
for(r = 1; r < matrix_size; r++) {
for(c = 0; c < col; c++) {
minor[r-1][c] = matrix[r][c];
}
for(c = col+1; c < matrix_size; c++) {
minor[r-1][c-1] = matrix[r][c];
}
}
// use "minor" matrix at this point to calculte
// its determinant
}

索引移位CCD_ 33是由于去除了第一行。

完全递归拉普拉斯展开

正如我之前提到的,行列式的拉普拉斯展开自然适用于递归算法。我将对您的设置进行一些更改,我不会使用堆栈分配的可变长度数组,而是使用堆分配的内存。由于如果不重用空间,则扩展具有指数空间需求,因此对于中等大小的矩阵,堆栈可能很快就会耗尽。因此,我将需要一个额外的参数来通过我称之为is_valid的intentit-out参数报告回内存分配失败。

您将认识到上面的矩阵复制过程有一些不同的名称和额外的指针取消引用,因为我在堆上使用n-by-n矩阵。我希望这不会导致太多的混乱。

#include <stdio.h>
#include <stdlib.h>

#define SQR(x) ((x)*(x))
int laplace_det(int matrix_size, const int (*mat)[][matrix_size], int *is_valid) {
// base cases
if(matrix_size == 1) 
return (*mat)[0][0];
if(matrix_size == 2)
return (*mat)[0][0] * (*mat)[1][1] - (*mat)[1][0] * (*mat)[0][1];
// recusive case, matrix_size > 2
// expansion indiscriminately along the first row
//
//   minor  matrix with i-th row and j-th column
//          removed for the purpose of calculating
//          the minor.
//   r, c   row and column index variables
//   col    current column in expansion
//   d      determinant accumulator
//
int r, c, col, d = 0;
int (*minor)[matrix_size-1][matrix_size-1] = calloc(SQR(matrix_size-1), sizeof(int));
if(!minor) {
*is_valid = 0;
return 0;
}
// expand along first row
for(col = 0; col < matrix_size; col++) {
// copy into minor matrix, left side then right side
for(r = 1; r < matrix_size; r++) {
for(c = 0; c < col; c++) {
(*minor)[r-1][c] = (*mat)[r][c];
}
for(c = col+1; c < matrix_size; c++) {
(*minor)[r-1][c-1] = (*mat)[r][c];
}
}
// calculate minor
int temp_d = laplace_det(matrix_size-1, minor, is_valid);
if(!is_valid) {
free(minor);
return 0;
}
d += (col & 1 ? -1 : 1) * (*mat)[0][col] * temp_d;
}
// free resources
free(minor);
return d;
}

示例驱动程序:

int main(void) {
int is_valid = 1;
int matrix[3][3] = {
{1, 2, 3},
{4, 5, 6},
{7, 8, 9}
};
int det_m = laplace_det(3, &matrix, &is_valid);
if(is_valid) {
printf("determinant %dn", det_m);
}
}

迭代方法

如果你想迭代地做同样的事情,你需要为所有较小尺寸的n-1子矩阵提供空间。正如递归情况所示,您可以对相同大小的所有子矩阵重用相同的空间,但不能对较小大小的矩阵使用该空间,因为每个矩阵必须生成较小的所有子数组,每个列一个。

由于矩阵的原始大小事先未知,因此很难实现以通用方式遍历所有这些矩阵,并且需要大量的记账,在递归的情况下,我们可以免费将这些值保留在它们各自的堆栈帧中。但我认为,将当前列选择器保留在大小为matrixSize的数组中,以及指向子矩阵的指针数组中,可以迭代重写。

我尝试使用递归方法来解决拉普拉斯展开。这能帮你吗

代码:

#include <stdio.h>

int determinant(int size,int det[][4])  // size & row of the square matrix
{
int temp[4][4],a=0,b=0,i,j,k;
int sum=0,sign;  /* sum will hold value of determinant of the current matrix */
if(size==2)
return (det[0][0]*det[1][1]-det[1][0]*det[0][1]);
sign=1;
for(i=0;i<size;i++)  // note that 'i' will indicate column no.
{
a=0;
b=0;
// copy into submatrix and recurse
for(j=1;j<size;j++) // should start from the next row !!
{
for(k=0;k<size;k++)
{
if(k==i) continue;
temp[a][b++]=det[j][k];
}
a++;
b=0;
}
sum+=sign*det[0][i]*determinant(size-1,temp);   // increnting row & decrementing size
sign*=-1;
}
return sum;
}
//Main function 
int main()
{
int i,j;
int det[4][4] = {{1, 0, 2, -1}, 
{3, 0, 0, 5}, 
{2, 1, 4, -3}, 
{1, 0, 5, 0} 
}; 
printf("%d",determinant(4,det));
}

干杯!

相关内容

  • 没有找到相关文章

最新更新