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?
One thought on “Matrix Inversion Using Swift and Accelerate”
Comments are closed.