使用LAPACK包装器估算Cython中LU分解的行列式 [英] Estimating determinant using LAPACK wrapper for LU decomposition in Cython

问题描述

我在这里定义了计算矩阵行列式的函数.但是有时候我得到了错误的信号.我通过此答案对函数进行了建模.

I define the function that calculates the determinant of a matrix here. But sometimes I get the wrong sign. I modeled my function from this answer.

from scipy.linalg.cython_lapack cimport dgetrf

cpdef double det_c(double[:, ::1] A, double[:, ::1] work, double[::1] ipiv):
    '''obtain determinant of float type square matrix A

    Notes
    -----
    As is, this function is not yet computing the sign of the determinant
    correctly, help!

    Parameters
    ----------
    A : memoryview (numpy array)
        n x n array to compute determinant of
    work : memoryview (numpy array)
        n x n array to use within function
    ipiv : memoryview (numpy array)
        length n vector use within function

    Returns
    -------
    detval : float
        determinant of matrix A
    '''

    cdef int n = A.shape[0], info
    work[...] = A

    dgetrf(&n, &n, &work[0,0], &n, &ipiv[0], &info)

    cdef double detval = 1.
    cdef int j

    for j in range(n):
        if j != ipiv[j]:
            detval = -detval*work[j, j]
        else:
            detval = detval*work[j, j]

    return detval

当我测试此功能并将其与np.linalg.det进行比较时，有时我得到了错误的符号.

When I test this function and compare it to np.linalg.det, sometimes I got the wrong sign.

>>> a = np.array([[1,2],[3,5.]])
>>> np.linalg.det(a)
>>> -1.0000000000000004
>>> det_c(a, np.zeros((2, 2)), np.zeros(2, dtype=np.int32))
>>> 1

其他时间，正确的标志.

Other times, the right sign.

>>> b = np.array([[1,2,3],[1,2,1],[5,6,1.]])
>>> np.linalg.det(b)
>>> -7.999999999999998
>>> det_c(a, np.zeros((3, 3)), np.zeros(3, dtype=np.int32))
>>> -8.0

推荐答案

dgetrf是Fortran子例程，并且Fortran使用基于1的索引，因此ipiv中的值在1到n之间(包括1和2). .为了解决这个问题，请在循环中从更改测试

dgetrf is a Fortran subroutine, and Fortran uses 1-based indexing, so the values in ipiv are between 1 and n (inclusive). To account for this, change the test in your loop from

        if j != ipiv[j]:

到

        if j != ipiv[j] - 1:

这篇关于使用LAPACK包装器估算Cython中LU分解的行列式的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！