Matrix Inversion Using Swift and Accelerate

I found out how to invert a matrix on SO, but I didn’t understand the solution, so I thought I’d talk more about it here. First of all, there isn’t a one-off inverse function in the Accelerate framework. You need to calculate the LU facotrization first using dgetrf_(), and then plug that data into dgetri_() to calculate the inverse.

The Solution

Assuming you’re only interested in the code, this is the code. This code assumes that the array m represents some matrix in row-major form, as in C or Python, not column-major form, as in Fortran. So m looks like,

| 1  2  3 |
| 0  1  4 |
| 5  6  0 | 

Furthermore, in the code, we’ll need to create and provide two additional arrays, pivots and workspace, and an integer error to the dgetrf_() and dgetri_() functions. We’ll pass all of the input to these functions as inout parameters. At the end of the day, we’ll return inMatrix as our inverted matrix.

import Accelerate

var m: [Double] = [ 1, 2, 3, 0, 1, 4, 5, 6, 0 ]

func invert(matrix : [Double]) -> [Double] {
    var inMatrix:[Double]       = matrix
    // Get the dimensions of the matrix. An NxN matrix has N^2
    // elements, so sqrt( N^2 ) will return N, the dimension
    var N:__CLPK_integer        = __CLPK_integer( sqrt( Double( matrix.count ) ) )
    // Initialize some arrays for the dgetrf_(), and dgetri_() functions
    var pivots:[__CLPK_integer] = [__CLPK_integer](count: Int(N), repeatedValue: 0)
    var workspace:[Double]      = [Double](count: Int(N), repeatedValue: 0.0)
    var error: __CLPK_integer   = 0
    // Perform LU factorization
    dgetrf_(&N, &N, &inMatrix, &N, &pivots, &error)
    // Calculate inverse from LU factorization
    dgetri_(&N, &inMatrix, &N, &pivots, &workspace, &N, &error)
    return inMatrix
}

var m_inv = invert( m )

print( m_inv )

This should give you something very close to,

"[-24.0, 18.0, 5.0, 20.0, -15.0, -4.0, -5.0, 4.0, 1.0]"

LU Factorization

According to the netlib source code for dgetrf(),

DGETRF computes an LU factorization of a general M-by-N matrix A using partial pivoting with row interchanges.

The factorization has the form

A = P * L * U

where P is a permutation matrix, L is lower triangular with unit diagonal elements (lower trapezoidal if m > n), and U is upper triangular (upper trapezoidal if m < n).

The arguments are given as follows,

Argument I/O Type Description
M input INTEGER The number of rows of the matrix A. M >= 0.
N input INTEGER The number of columns of the matrix A. N >= 0
A input/output DOUBLE PRECISION array, dimension (LDA,N) On entry, the M-by-N matrix to be factored. On exit, the factors L and U from the factorization A = P*L*U; the unit diagonal elements of L are not stored.
LDA input INTEGER The leading dimension of the array A. LDA >= max(1,M).
IPIV output INTEGER array, dimension (min(M,N)) The pivot indices; for 1 <= i <= min(M,N), row i of the matrix was interchanged with row IPIV(i).
INFO output INTEGER = 0: successful exit
< 0: the i-th argument had an illegal value
> 0: U(i,i) is exactly zero. The factorization has been completed, but the factor U is exactly singular, and division by zero will occur if it is used to solve a system of equations.

Matrix Inversion

According to the netlib source code for dgetri()

DGETRI computes the inverse of a matrix using the LU factorization computed by DGETRF.

This method inverts U and then computes inv(A) by solving the system inv(A)*L = inv(U) for inv(A).

The arguments are given as follows,

Argument I/O Type Description
N input INTEGER The order of the matrix A. N >= 0.
A input/output DOUBLE PRECISION array, dimension (LDA,N) On entry, the factors L and U from the factorization A = P*L*U as computed by DGETRF. On exit, if INFO = 0, the inverse of the original matrix A.
LDA input INTEGER The leading dimension of the array A. LDA >= max(1,N).
IPIV input INTEGER array, dimension (N) The pivot indices from DGETRF; for 1<=i<=N, row i of the matrix was interchanged with row IPIV(i).
WORK input/output DOUBLE PRECISION array, dimension (MAX(1,LWORK)) On exit, if INFO=0, then WORK(1) returns the optimal LWORK.
LWORK input INTEGER The dimension of the array WORK. LWORK >= max(1,N). For optimal performance LWORK >= N*NB, where NB is the optimal blocksize returned by ILAENV.
If LWORK = -1, then a workspace query is assumed; the routine only calculates the optimal size of the WORK array, returns this value as the first entry of the WORK array, and no error message related to LWORK is issued by XERBLA.
INFO output INTEGER = 0: successful exit
< 0: the i-th argument had an illegal value
> 0: U(i,i) is exactly zero; the matrix is singular and its inverse could not be computed.

References

Matrix Inversion in Swift using Accelerate Framework

How to perform matrix inverse operation using the accelerate framework?

Using Accelerate (CLAPACK) to solve an augmented matrix?

dgetrf()

dgetri()

One thought on “Matrix Inversion Using Swift and Accelerate”

Leave a Reply

Your email address will not be published. Required fields are marked *