Как написать функции для вычисления обратной матрицы и решить в Cython, используя Lapack, без объектов Python? - PullRequest
Новая
       40

Как написать функции для вычисления обратной матрицы и решить в Cython, используя Lapack, без объектов Python?

1 голос
/ 09 марта 2019

Я пытаюсь написать функции Cython для вычисления инверсий и получения решений для линейных систем.Я знаю, что это то, что делает Numpy.Но мне нужно использовать эти функции в другом цикле, который я пишу на Cython.Итак, я не хочу использовать объекты Python в цикле.

Я написал обратную функцию и функцию решения, и они иногда работают.Но в других случаях они терпят неудачу.Я думаю, что я могу неправильно понять что-то фундаментальное.Во всяком случае, здесь MWE.Все файлы находятся в одном каталоге.Сначала файл mopp.pxd: 

#cython: boundscheck=False, wraparound=False, nonecheck=False
from scipy.linalg.cython_lapack cimport dgetri, dgetrf, dgesv
cimport cython

cdef inline void solve_c(double[:, ::1] A, double[:, ::1] B, double[:, ::1] out,
    int[::1] ipiv, double[:, ::1] Awork):
    '''solve system A*X=B [WARNING: Not working for B with 2+ columns]

    Parameters
    ----------
    A : memoryview
        n x n
    B : memoryview
        n x r
    out : memoryview
        n x r, for output storage
    ipiv : memoryview
        length n integer vector
    Awork : memoryview
        n x n vector to use for work
    '''

    cdef int LDA = A.shape[0], N = A.shape[1]
    cdef int LDB = B.shape[0], NRHS = B.shape[1]
    cdef int info

    Awork[...] = A
    out[...] = B

    dgesv(&N, &NRHS, &Awork[0,0], &LDA, &ipiv[0], &out[0,0], &LDB, &info)

cdef inline void inv_c(double[:, ::1] A, double[:, ::1] B, 
                          double[:, ::1] work, int[::1] ipiv):
    '''invert float type square matrix A [WARNING: Not working randomly?]

    Parameters
    ----------
    A : memoryview (numpy array)
        n x n array to invert
    B : memoryview (numpy array)
        n x n array to use within the function, function
        will modify this matrix in place to become the inverse of A
    work : memoryview (numpy array)
        n x n array to use within the function
    ipiv : memoryview (numpy array)
        length n array to use within function
    '''

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

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

    dgetri(&n, &B[0,0], &n, &ipiv[0], &work[0,0], &lwork, &info)

Затем файл, содержащий обертки, чтобы сделать эти функции доступными для тестирования Python, под названием mopp_wrappers.pyx: 

#cython: boundscheck=False, wraparound=False, nonecheck=False, cdivision=True
from mopp cimport inv_c, solve_c
import numpy as np

def inv(A, B, work, ipiv):

    inv_c(A, B, work, ipiv)

def solve(A, B, out, ipiv, Awork):

    solve_c(A, B, out, ipiv, Awork)

И у меня естьсвязанный файл установки setup.py: 

from distutils.core import setup
from Cython.Build import cythonize

setup(name="mopp_wrappers", ext_modules=cythonize('mopp_wrappers.pyx', annotate=True))

Затем я компилирую, запустив в терминале следующее: 

python setup.py build_ext --inplace

Наконец, я проверяю функции в Python с test.py.Я просто не понимаю, почему функции работают в некоторых сценариях, но не в других. 

import numpy as np
from mopp_wrappers import inv, solve

def mysolve(A, b):
    out = np.zeros((A.shape[0], b.shape[1]))
    ipiv = np.zeros(3, dtype=np.int32)
    Awork = np.zeros(A.shape)

    solve(A, b, out, ipiv, Awork)

    return out

def myinv(A):
    B = np.zeros(A.shape)
    work = np.zeros(A.shape)
    ipiv = np.zeros(A.shape[0], dtype=np.int32)

    inv(A, B, work, ipiv)

    return B

np.random.seed(100)

M = np.random.multivariate_normal(mean=np.zeros(3), cov=np.eye(3), size=20)

A = M.T.dot(M)

# inverse is the same
print('Successful computation of inverse:')
print(np.linalg.inv(A))
print(myinv(A))
print('')

# inverse is different
B = np.zeros(A.shape)
work = np.zeros(A.shape)
ipiv = np.zeros(A.shape[0], dtype=np.int32)

inv(A, B, work, ipiv)

print('Failure to compute inverse:')
print(B)
print(np.linalg.inv(A))
print('')


# solution is the same
b = np.array([1,2,3]).reshape((-1, 1)).astype(float)

print('Successful computation of solution, 1-dim b matrix')
print(mysolve(A, b))
print(np.linalg.solve(A, b))
print('')

# solution is different
b = np.hstack([b, b])
print('Failure to compute solution, 2-dim b matrix')
print(mysolve(A, b))
print(np.linalg.solve(A, b))

, где я получаю следующий вывод: 

Successful computation of inverse:
[[ 0.0587434   0.01003223  0.01867754]
 [ 0.01003223  0.06387784 -0.00035565]
 [ 0.01867754 -0.00035565  0.06977158]]
[[ 0.0587434   0.01003223  0.01867754]
 [ 0.01003223  0.06387784 -0.00035565]
 [ 0.01867754 -0.00035565  0.06977158]]

Failure to compute inverse:
[[ 1.91799368e+01 -1.58548362e-01 -2.68503736e-01]
 [-3.04094756e+00  1.56553248e+01  5.09740667e-03]
 [-5.14988469e+00  7.98015569e-02  1.43324831e+01]]
[[ 0.0587434   0.01003223  0.01867754]
 [ 0.01003223  0.06387784 -0.00035565]
 [ 0.01867754 -0.00035565  0.06977158]]

Successful computation of solution, 1-dim b matrix
[[0.13484049]
 [0.13672096]
 [0.22728098]]
[[0.13484049]
 [0.13672096]
 [0.22728098]]

Failure to compute solution, 2-dim b matrix
[[0.10613072 0.07319877]
 [0.15786505 0.20361612]
 [0.21063103 0.24560286]]
[[0.13484049 0.13484049]
 [0.13672096 0.13672096]
 [0.22728098 0.22728098]]
...