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

How do I implement the Laplace expansion algorithm in 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。

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

我试图用递归方法解决拉普拉斯展开。希望对您有所帮助

代码:

#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));
}

干杯!

Those values have to be multiplied by the matrix that it's left when I remove the row and column where that value is.

您必须乘以余数矩阵,其条目是删除第 i 行和第 j 列时留下的矩阵的行列式

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

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

你的方法

让我们看看您提供的这段代码。它使我无法理解您要实现的目标。以迭代 col:

的最内层循环中的语句为例
matrixFirstRowValues[row + (matrixSize - i)] = matrix[1][col + (matrixSize - i)];

只要最内层循环中的 col 发生变化,rowi 都是固定的,因此您正在从(显然)第二行分配和重新分配 matrixmatrixFirstRowValues 中的相同条目。不仅如此,您从超出列范围的索引范围 (matrixSize-i) .. (2*matrixSize - (i+1)) 分配,除非 i == matrixSize,这仅适用于 i.

的第一个值

正如我之前提到的,最后你得到的不仅仅是一个 2×2 矩阵,而是 n!/2

复制第i行第j列除外

查看删除了第 i 行和第 j 列的矩阵,您最终得到四个子矩阵(其中一些可能为空)。限制自己沿着第一行展开,然后你只处理两个子矩阵(其中一些可能是空的)。您可以使用两个循环,一个用于第 j 列左侧和右侧的矩阵 - 或者,如上一个答案中所建议的那样 - 选择跳过第 j 列使用continue 循环循环而不更新目标列索引。如果 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
}

索引偏移 r-1 是由于删除了第一行。

一个完整的递归拉普拉斯展开

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

的 intent out 参数报告内存分配失败。

你会发现上面的矩阵复制过程有一点不同的名字和一个额外的指针取消引用,因为我在堆上操作 n×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 %d\n", det_m);
    }
}

迭代方法

如果您想迭代地做同样的事情,您将需要为所有 n-1 个更小的子矩阵提供 space。正如递归案例所示,您可以对 相同 大小的所有子矩阵重用相同的 space,但不能将 space 用于较小大小的矩阵,因为每个矩阵必须产生所有尺寸小一的子矩阵,每列一个。

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