Copyright (c) 2020 Steven Obua
License: MIT License
LANumerics is a Swift package for doing numerical linear algebra.
The package depends on Swift Numerics, as it supports both real and complex numerics for both Float
and Double
precision in a uniform way.
Under the hood it relies on the Accelerate
framework for most of its functionality, in particular BLAS
and LAPACK
, and also vDSP
.
Examining the current tests provides a good starting point beyond this README.
LANumerics is a normal Swift package and can be added to your app in the usual way.
After adding it to your app, import LANumerics
(and also Numerics
if you use complex numbers).
You can try out if everything works fine by running
import LANumerics
let A : Matrix<Float> = Matrix(columns: [[1, 2, 3, 4], [5, 6, 7, 8], [9, 10, 11, 12]])
print("A: \(A)")
which should output something like
A: 4x3-matrix:
β1.0 5.0 9.0 β
β2.0 6.0 10.0β
β3.0 7.0 11.0β
β4.0 8.0 12.0β
The LANumeric
protocol denotes the type of numbers on which LANumerics operates. It is implemented by the following types:
Float
Double
Complex<Float>
Complex<Double>
Most functionality of LANumerics is generic in LANumeric
, e.g. constructing matrices and computing with them, solving a system of linear equations, or computing the singular value decomposition of a matrix.
The main work horse of LANumerics is the Matrix
type. For convenience there is also a Vector
type, but this is just a typealias for normal Swift arrays.
The expression Matrix([1,2,3])
constructs the matrix:
3x1-matrix:
β1.0β
β2.0β
β3.0β
The expression Matrix(row: [1, 2, 3])
constructs the matrix:
1x3-matrix:
(1.0 2.0 3.0)
The expression Matrix<Float>(rows: 2, columns: 3)
constructs a matrix consisting only of zeros:
2x3-matrix:
β0.0 0.0 0.0β
β0.0 0.0 0.0β
The expression Matrix(repeating: 1, rows: 2, columns: 3)
constructs a matrix consisting only of ones:
2x3-matrix:
β1.0 1.0 1.0β
β1.0 1.0 1.0β
Given the two vectors v1
and v2
let v1 : Vector<Float> = [1, 2, 3]
let v2 : Vector<Float> = [4, 5, 6]
we can create a matrix from columns Matrix(columns: [v1, v2])
:
3x2-matrix:
β1.0 4.0β
β2.0 5.0β
β3.0 6.0β
or rows Matrix(rows: [v1, v2])
:
2x3-matrix:
β1.0 2.0 3.0β
β4.0 5.0 6.0β
It is also legal to create matrices with zero columns and/or rows, like Matrix(rows: 2, columns: 0)
or Matrix(rows: 0, columns: 0)
.
Swift supports simd
vector and matrix operations. LANumerics plays nice with simd
by providing conversion functions to and from simd
vectors and matrices. For example, starting from
import simd
import LANumerics
let m = Matrix(rows: [[1, 2, 3], [4, 5, 6]])
print("m: \(m)")
with output
m: 2x3-matrix:
β1.0 2.0 3.0β
β4.0 5.0 6.0β
we can convert m
into a simd
matrix s
via
let s = m.simd3x2
print(s)
resulting in the output
simd_double3x2(columns: (SIMD2<Double>(1.0, 4.0), SIMD2<Double>(2.0, 5.0), SIMD2<Double>(3.0, 6.0)))
Note that simd
reverses the role of row and column indices compared to LANumerics
(and usual mathematical convention).
We can also convert s
back:
print(Matrix(s) == m)
will yield the output true
.
Matrix elements and submatrices can be accessed using familiar notation. Given
import simd
import LANumerics
var m = Matrix(rows: [[1, 2, 3], [4, 5, 6], [7, 8, 9]])
print(m)
with output
3x3-matrix:
β1.0 2.0 3.0β
β4.0 5.0 6.0β
β7.0 8.0 9.0β
we can access the element at row 2 and column 1 via
m[2, 1]
which yields 8.0
. We can also set the element at row 2 and column 1 to some value:
m[2, 1] = 0
print(m)
The output of running this is
3x3-matrix:
β1.0 2.0 3.0β
β4.0 5.0 6.0β
β7.0 0.0 9.0β
We can also access submatrices of m
, for example its top-left and bottom-right 2x2 submatrices:
print(m[0 ... 1, 0 ... 1])
print(m[1 ... 2, 1 ... 2])
This will print
2x2-matrix:
β1.0 2.0β
β4.0 5.0β
2x2-matrix:
β5.0 6.0β
β0.0 9.0β
Finally, using the same notation, we can overwrite submatrices of m
:
m[0 ... 1, 0 ... 1] = m[1 ... 2, 1 ... 2]
print(m)
This overwrites the top-left 2x2 submatrix of m
with its bottom-right 2x2 submatrix, yielding:
3x3-matrix:
β5.0 6.0 3.0β
β0.0 9.0 6.0β
β7.0 0.0 9.0β
LANumerics supports common operations on matrices and vectors, among them:
transpose
and adjoint
In the following, assume the context
import Numerics
import LANumerics
let u = Matrix<Complex<Float>>(rows: [[1, 2 * .i], [3, 4 * .i + 1]])
let v = Matrix<Complex<Float>>(rows: [[.i, 0], [0, 1 + 1 * .i]])
print("u : \(u)\n")
print("v : \(v)\n")
which has output
u : 2x2-matrix:
β1.0 2.0i β
β3.0 1.0 + 4.0iβ
v : 2x2-matrix:
β1.0i 0.0 β
β0.0 1.0 + 1.0iβ
For real matrices, transpose
and adjoint
have the same meaning, but for complex matrices the adjoint
is the element-wise conjugate of the transpose
. Executing
print("u.transpose : \(u.transpose)\n")
print("u.adjoint : \(u.adjoint)\n")
thus yields
u.transpose : 2x2-matrix:
β1.0 3.0 β
β2.0i 1.0 + 4.0iβ
u.adjoint : 2x2-matrix:
β1.0 3.0 β
β-2.0i 1.0 - 4.0iβ
The adjoint
has the advantage over the transpose
that many properties involving the adjoint generalize naturally from real matrices to complex matrices. Therefore there is the shortcut notation
uβ²
for u.adjoint
.
Note that β²
is the unicode character "Prime" U+2032
. You can use for example Ukelele to make the input of that character smooth. Other alternatives are configuring the touchbar of your macbook, or using a configurable keyboard like Stream Deck.
Multiplying u
and v
is done via the expression u * v
. Running print("u * v: \(u * v)")
results in
u * v: 2x2-matrix:
β1.0i -2.0 + 2.0iβ
β3.0i -3.0 + 5.0iβ
Instead of uβ² * v
one can also use the equivalent, but faster expression u β²* v
:
print("uβ² * v: \(uβ² * v)\n")
print("u β²* v: \(u β²* v)\n")
yields
uβ² * v: 2x2-matrix:
β1.0i 3.0 + 3.0iβ
β2.0 5.0 - 3.0iβ
u β²* v: 2x2-matrix:
β1.0i 3.0 + 3.0iβ
β2.0 5.0 - 3.0iβ
Similarly, it is better to use u *β² v
than u * vβ²
, and u β²*β² v
instead of uβ² * vβ²
.
We will view u
and v
as vectors u.vector
and v.vector
now, where .vector
corresponds to a column-major order of the matrix elements:
u.vector: [1.0, 3.0, 2.0i, 1.0 + 4.0i]
v.vector: [1.0i, 0.0, 0.0, 1.0 + 1.0i]
(Actually, in the above, we used u.vector.toString()
and v.vector.toString()
for better formatting of complex numbers. We will also do so below where appropriate without further mentioning it.)
The dot product of u.vector
and v.vector
results in
u.vector * v.vector: -3.0 + 6.0i
Another vector product is
u.vector β²* v.vector: 5.0 - 2.0i
which corresponds to
u.vectorβ² * v.vector: [5.0 - 2.0i]
Furthermore, there is
u.vector *β² v.vector: 4x4-matrix:
β-1.0i 0.0 0.0 1.0 - 1.0iβ
β-3.0i 0.0 0.0 3.0 - 3.0iβ
β2.0 0.0 0.0 2.0 + 2.0iβ
β4.0 - 1.0i 0.0 0.0 5.0 + 3.0iβ
which is equivalent to
Matrix(u.vector) * v.vectorβ²: 4x4-matrix:
β-1.0i 0.0 0.0 1.0 - 1.0iβ
β-3.0i 0.0 0.0 3.0 - 3.0iβ
β2.0 0.0 0.0 2.0 + 2.0iβ
β4.0 - 1.0i 0.0 0.0 5.0 + 3.0iβ
Element-wise operations like .+
, .-
, .*
and ./
are supported on both vectors and matrices, for example:
u .* v : 2x2-matrix:
β1.0i 0.0 β
β0.0 -3.0 + 5.0iβ
The Matrix
type supports functional operations like map
, reduce
and combine
. These come in handy when performance is not that important, and there is no accelerated equivalent available (yet?).
For example, the expression
u.reduce(0) { x, y in max(x, y.magnitude) }
results in the value 4
. In this case it is better though to use the equivalent expression u.infNorm
instead.
You can solve a system of linear equations like
by converting it into matrix form A * u = b
and solving it for u
:
import LANumerics
let A = Matrix<Double>(rows: [[7, 5, -3], [3, -5, 2], [5, 3, -7]])
let b : Vector<Double> = [16, -8, 0]
let u = A.solve(b)!
print("A: \(A)\n")
print("b: \(b)\n")
print("u: \(u.toString(precision: 1))")
This results in the output:
A: 3x3-matrix:
β7.0 5.0 -3.0β
β3.0 -5.0 2.0 β
β5.0 3.0 -7.0β
b: [16.0, -8.0, 0.0]
u: [1.0, 3.0, 2.0]
Therefore the solution is x=1, y=3, z=2.
This example is actually plugged from an article which describes how to use the Accelerate
framework directly and without a nice library like LANumerics
π.
You can solve for multiple right-hand sides simultaneously. For example, you can compute the inverse of A
like so:
let Id : Matrix<Double> = .eye(3)
print("Id: \(Id)\n")
let U = A.solve(Id)!
print("U: \(U)\n")
print("A * U: \(A * U)")
This results in the output
Id: 3x3-matrix:
β1.0 0.0 0.0β
β0.0 1.0 0.0β
β0.0 0.0 1.0β
U: 3x3-matrix:
β0.11328125 0.1015625 -0.019531249999999997β
β0.12109375 -0.13281250000000003 -0.08984375 β
β0.1328125 0.01562499999999999 -0.1953125 β
A * U: 3x3-matrix:
β1.0 -1.0755285551056204e-16 0.0β
β0.0 1.0 0.0β
β0.0 -4.163336342344337e-17 1.0β
Extending the above example, you can also solve it using least squares approximation instead:
print(A.solveLeastSquares(b)!.toString(precision: 1))
results in the same solution [1.0, 3.0, 2.0]
.
Least squares is more general than solving linear equations directly, as it can also deal with situations where you have more equations than variables, or less equations than variables. In other words, it can also handle:
A
There is a shorthand notation available for the expression A.solveLeastSquares(b)!
:
A β b
This also works for simultaneously solving for multiple right-hand sides as before, the inverse of A
can therefore also be computed using the expression A β .eye(3)
:
A β .eye(3): 3x3-matrix:
β0.11328124999999997 0.10156250000000001 -0.01953124999999994β
β0.12109375 -0.13281250000000006 -0.08984374999999999β
β0.13281249999999997 0.015624999999999993 -0.19531249999999994β
Note that β
is the unicode character "Set Minus" U+2216
. The same advice for smooth input of this character applies as for the input of β²
earlier.
In addition, there is the operator β²β
which combines taking the adjoint and solving via least squares. Therefore, to compute the inverse of Aβ²
, you could better write A β²β .eye(3)
instead of Aβ² β .eye(3)
:
A β²β .eye(3): 3x3-matrix:
β0.11328124999999996 0.12109375 0.13281249999999994 β
β0.10156250000000001 -0.13281250000000003 0.01562499999999999 β
β-0.01953124999999994 -0.08984374999999999 -0.19531249999999994β
The inverse of A
can more concisely also be obtained via A.inverse!
.
The following matrix decompositions are currently supported:
A
:
A.svd()
A
(that is for real or complex matrices for which A == Aβ²
):
A.eigen()
A
:
A.schur()
link |
Stars: 35 |
Last commit: 2 years ago |
Swiftpack is being maintained by Petr Pavlik | @ptrpavlik | @swiftpackco | API | Analytics