Linear algebra library based on LAPACK
Linear algebra library in Crystal, uses LAPACK.
Killing SciPy, one module at a time.
Install LAPACK and BLAS. sudo apt install libopenblas-base liblapack3
for Ubuntu, sudo pacman -S lapack
(for better performance use openblas-lapack
package from AUR) for Arch. For Windows you need libopenblas.dll (and libopenblas.lib) from https://github.com/xianyi/OpenBLAS/releases
(for Ubuntu 18) it seems package doesn't create symlink, so use
sudo ln -s /usr/lib/lapack/liblapack.so.3 /usr/lib/liblapack.so
sudo ln -s /usr/lib/openblas-base/libblas.so.3 /usr/lib/libcblas.so
shard.yml
:dependencies:
linalg:
github: konovod/linalg
shards install
require "linalg"
Basic type aliases are
Complex consisting of two Float32 isn't supported for now (it is easy, but I'm not sure if it's useful).
Types with prefix G (GMat, GMat32, GMatComplex) are for actually allocated matrices, others are automatically converted to them when needed.
#suggested to don't prefix LA:: everywhere
include LA
# create matrix from array of arrays (or tuple... everything Indexable)
m = GMat[
[1, 2, 3],
[4, 5, 6],
[7, 8, 9],
[10, 11, 12],
]
# or using block
m = GMat32.new(3, 4) { |i, j| i*3 + j + 1 }
# or using one of other ways, check "spec" directory
# do basic arithmetics
a = Mat.eye(3)
pp 2 * a - Mat.diag([2, 2, 2]) == Mat.zeros(3, 3) # => true
# basic algebra
a = Mat.rand(5, 5) + 2 * Mat.identity(5)
pp (a.inv * a - Mat.identity(5)).norm < 1e-6
b = Mat.rand(5, 1)
x = LA.solve(a, b) # or a.solve(b)
pp (a*x - b).norm < 1e-6
m = GMat[[-2, 4, 1], [2, -4, 1], [1, 1, 1]]
pp m.eigvals # => [-6.0, -1.0, 2.0]
# extract submatrices (memory isn't copied as they reference to basic matrix)
m = GMat[
[1, 2, 3],
[4, 5, 6],
[7, 8, 9],
]
pp m.columns[2] # LA::SubMatrix(Float64) (3x1, None):
# [3.0]
# [6.0]
# [9.0]
x = m[1..1, 1..2]
pp x # => [5.0, 6.0]
x[0, 0] = 0 # m[1,1] is now 0 (questionable feature? maybe should be ##[]! for modifiable submatrices and ##[] for CoW?)
y = x.clone # now y is a separate matrix
y[0, 0] = 1 # m[1,1] is still 0
pp m[1, 1]
other present features:
Mat#svd
or Mat#svdvals
for just values)Mat#lu
)# to just get P L U matrices
p, l, u = a.lu
# to get them in compact form and use for solving linear equations:
a = GMat32[
[2, 4],
[2, 8]
]
lu = a.lu_factor # lu is LUMatrix(T) - immutable object that can return it's content and solve systems
puts lu.solve(GMat32[[2], [4]])
LA.solvels
to just get decision or LA.lstsq
to also get rank and singular values (TODO - and residues))#cholesky
, #cholesky!
, #cho_solve
)hessenberg
formqr
, rq
, lq
, ql
decompositionsschur
and qz
(generalized schur) decompositioneigs(a, b, ...)
)pascal
or toeplitz
(check scipy.md for a full list)There is also concept of Mat#flags
that represent properties of matrix (symmetric, positive definite etc), they are used to automatically select faster algorithms from LAPACK. Flags are partially enforced by runtime checks, with the possibility of user override. For example, if we say that a.assume!(MatrixFlags::Symmetric)
then a.transpose
or a + Mat.diag(*a.size)
will also have this flag, so the LAPACK routines for symmetrical matrices will be used. In fact, a.transpose
will return matrix clone as for symmetric matrices A=A'.
Supported flags:
enum MatrixFlags
Symmetric
Hermitian
PositiveDefinite
Orthogonal
UpperTriangular
LowerTriangular
Triangular = UpperTriangular | LowerTriangular
NOTE for complex matrices Orthogonal
flag means Unitary
.
Main functions for flags are:
a.assume!(flag) # sets matrix flag without check, can lead to incorrect results if matrix do not have corresponding property.
a.detect?(flag) # checks if matrix has property, if yes sets the flag. Returns true if check positive
a.detect(flag) # same as `detect?`, but returns matrix a
a.detect # detect all possible flags
a.flags # returns matrix flags
Most operations - matrix addition, multiplication, inversion, transposition and decompositions correctly update flags, but any direct access like a[i,j] = 0
or a.map!{|v| v+1}
resets flags to None
, so use a.detect
after them if you need to preserve flags (or a.assume!(f)
if detection is too slow).
saving/loading from files
ways to evade allocations during calculations
Matrix exponent and trigonometric
other matrix functions (fractional power and logarithm)
Banded matrices
Column-major storage (optional?)
Other missing features from LAPACK (selectable and orderable eigenvalues, generalized decompositions, new routines etc)
Sparse matrices
Other missing features from scipy.linalg (lyapunov/ricatti/sylvester equations, other things i don't know algorithms for)
Support multithreading\reenterability
Allow matrices of unsupported types (use check_type
only when actually needed)
T.new(value)
, maybe this requirement can be relaxedmap
signature has to be changed?arange
, open-ended ranges in submatricesGMat.max
fix