PyArmadillo: a streamlined linear algebra library for Python

PyArmadillo is a linear algebra library for the Python language, with the aim of closely mirroring the programming interface of the widely used Armadillo C++ library, which in turn is deliberately similar to Matlab. PyArmadillo hence facilitates algorithm prototyping with Matlab-like syntax directly in Python, and relatively straightforward conversion of PyArmadillo-based Python code into performant Armadillo-based C++ code. The converted code can be used for purposes such as speeding up Python-based programs in conjunction with pybind11, or the integration of algorithms originally prototyped in Python into larger C++ codebases. PyArmadillo provides objects for matrices and cubes, as well as over 200 associated functions for manipulating data stored in the objects. Integer, floating point and complex numbers are supported. Various matrix factorisations are provided through integration with LAPACK, or one of its high performance drop-in replacements such as Intel MKL or OpenBLAS. PyArmadillo is open-source software, distributed under the Apache 2.0 license; it can be obtained at https://pyarma.sourceforge.io or via the Python Package Index in precompiled form.

PyArmadillo is aimed at: (i) users that prefer compact Matlab-like syntax rather than the somewhat more verbose syntax provided by NumPy/SciPy [6,17], and (ii) users that would like a straightforward conversion path to performant C++ code. More specifically, PyArmadillo aims to closely mirror the programming interface of the Armadillo library, thereby facilitating the prototyping of algorithms with Matlab-like syntax directly in Python. Furthermore, PyArmadillo-based Python code can be easily converted into high-performance Armadillo-based C++ code. Due to the similarity of the programming interfaces, the risk of introducing bugs in the conversion process is considerably reduced. Moreover, conversion into C++ based code allows taking advantage of expression optimisation performed at compile-time by Armadillo, resulting in further speedups. The resulting code can be used in larger C++ programs, or used as a replacement of performance critical parts within a Python program with the aid of the pybind11 interface layer [9].

Functionality
PyArmadillo provides matrix objects for several distinct element types: integers, single-and double-precision floating point numbers, as well as complex numbers. In addition to matrices, PyArmadillo also has support for cubes (3 dimensional arrays), where each cube can be treated as an ordered set of matrices. Over 200 functions are provided for manipulating data stored in the objects, covering the following areas: fundamental arithmetic operations, contiguous and non-contiguous submatrix views, diagonal views, element-wise functions, scalar/vector/matrix valued functions of matrices, statistics, signal processing, storage of matrices in files, matrix decompositions/factorisations, matrix inverses, and equation solvers. PyArmadillo matrices and cubes are convertible to/from NumPy arrays, allowing users to tap into the wider Python data science ecosystem, including plotting tools such as Matplotlib [7].
An overview of the available functionality in PyArmadillo (as of version 0.500) is given in Tables 1 through to 8.  Table 1 briefly describes the member functions and variables of the main matrix class; Table 2 lists the main subset of overloaded Python operators; Table 3 outlines matrix decompositions and equation solvers;  Table 4 overviews functions for generating vectors with various sequences; Table 5 lists the main forms of general functions of matrices; Table 6 lists various element-wise functions of matrices; Table 7 summarises the set of functions and classes focused on statistics; Table 8 lists functions for signal and image processing. Lastly, Table 9 provides examples of Matlab syntax and conceptually corresponding PyArmadillo syntax.  Tables 1 to 8 for an overview of the available functionality.

Implementation & License
PyArmadillo relies on pybind11 [9] for interfacing C++ and Python, as well as on Armadillo [13,14] for the underlying C++ implementation of matrix objects and associated functions. Due to its expressiveness and relatively straightforward use, pybind11 was selected over other interfacing approaches such as Boost.Python [1] and manually writing C++ extensions for Python. In turn, Armadillo interfaces with low-level routines in BLAS and LAPACK [2], where BLAS is used for matrix multiplication, and LAPACK is used for various matrix decompositions/factorisations and equation solvers. As the low-level routines in BLAS and LAPACK are considered as a de facto standard for numerical linear algebra, it is possible to use high performance drop-in replacements such as Intel MKL [8] and OpenBLAS [18]. PyArmadillo is open-source software available under the permissive Apache License 2.0 [3], making it useful in both open-source and proprietary (closed-source) contexts [16]. access the i-th element, assuming a column-by-column layout [r, c] access the element at row r and column c .in range(i) test whether the i-th element can be accessed .in range(r, c) test whether the element at row r and column c can be accessed .reset() set the number of elements to zero .copy size(A) set the size to be the same as matrix A .set size(n rows, n cols) change size to specified dimensions, without preserving data (fast) .reshape(n rows, n cols) change size to specified dimensions, with elements copied column-wise (slow) .resize(n rows, n cols) change size to specified dimensions, while preserving elements & their layout (slow) .ones(n rows, n cols) set all elements to one, optionally first resizing to specified dimensions .zeros(n rows, n cols) as above, but set all elements to zero .randu(n rows, n cols) as above, but set elements to uniformly distributed random values in [0,1] interval .randn(n rows, n cols) as above, but use a Gaussian/normal distribution with µ = 0 and σ = 1 .fill(k) set all elements to be equal to k .for each(lambda val : ...) for each element, pass its value to a lambda function .is empty() test whether there are no elements read/write access to matrix elements corresponding to the specified indices .swap rows(p, q) swap the contents of specified rows .swap cols(p, q) swap the contents of specified columns .insert rows(row, X) insert a copy of X at the specified row .insert cols(col, X) insert a copy of X at the specified column .shed rows(first row, last row) remove the specified range of rows .shed cols(first col, last col) remove the specified range of columns .min() return minimum value .max() return maximum value .index min() return index of minimum value .index max() return index of maximum value  Cholesky decomposition of symmetric positive-definite matrix X eig sym(X) eigen decomposition of a symmetric/hermitian matrix X eig gen(X) eigen decomposition of a general (non-symmetric/non-hermitian) square matrix X eig pair(A, B) eigen decomposition for pair of general square matrices A and B inv(X) inverse of a square matrix X inv sympd(X) inverse of symmetric positive definite matrix X lu(L, U, P, X) lower-upper decomposition of X, such that PX = LU and X = P'LU null(X) orthonormal basis of the null space of matrix X orth(X) orthonormal basis of the range space of matrix X pinv(X) Moore-Penrose pseudo-inverse of a non-square matrix X qr(Q, R, X) QR decomposition of X, such that QR = X qr econ(Q, R, X) economical QR decomposition qz(AA, BB, Q, Z, A, B) generalised Schur decomposition for pair of general square matrices A and B schur(X) Schur decomposition of square matrix X solve(A, B) solve a system of linear equations AX = B, where X is unknown svd(X) singular value decomposition of X svd econ(X) economical singular value decomposition of X syl(X) Sylvester equation solver  vector with m elements, with random integers sampled without replacement from 0 to n-1 randu(n rows, n cols) matrix with random values from a uniform distribution in the [0,1] interval randn(n rows, n cols) matrix with random values from a normal/Gaussian distribution with µ = 0 and σ = 1 zeros(n rows, n cols) matrix with all elements set to zero ones(n rows, n cols) matrix with all elements set to one  p-norm of matrix A, with p = 1, 2, · · ·, or p = "-inf", "inf", "fro" normalise(A, p, dim) return the normalised version of A, with each column or row normalised to unit p-norm prod(A, dim) product of elements in each column or row of matrix A rank(A) rank of matrix A rcond(A) estimate the reciprocal of the condition number of square matrix A real(C) extract the real part of complex matrix C repmat (A, p, q) replicate matrix A in a block-like fashion, resulting in p by q blocks of matrix A reshape(A, r, c) create matrix with r rows and c columns by copying elements from A column-wise resize(A, r, c) create matrix with r rows and c columns by copying elements and their layout from A shift(A, n, dim) copy matrix A with the elements shifted by n positions in each column or row shuffle(A, dim) copy matrix A with elements shuffled in each column or row size(A) obtain the dimensions of matrix A sort(A, direction, dim) copy A with elements sorted (in ascending or descending direction) in each column or row sort index(A, direction) generate a vector of indices describing the sorted order of the elements in matrix A sqrtmat(A) complex square root of square matrix A sum(A, dim) sum of elements in each column or row of matrix A sub2ind(size(A), row, col) convert subscript notation (row,col) to a linear index, using the size of matrix  Table 6: Element-wise functions: matrix B is produced by applying a function to each element of matrix A.
trignometric function, where trig is one of: cos, acos, cosh, acosh, sin, asin, sinh, asinh, tan, atan, tanh, atanh generate matrix of histogram counts for each column or row of A, using given bin centers histc(A, edges, dim) generate matrix of histogram counts for each column or row of A, using given bin edges kmeans(means, A, k, ...) cluster column vectors in matrix A into k disjoint sets, storing the set centers in means princomp(A) principal component analysis of matrix A running stat class for running statistics of a continuously sampled one dimensional signal running stat vec class for running statistics of a continuously sampled multi-dimensional signal mean(A, dim) find the mean in each column or row of matrix A median(A, dim) find the median in each column or row of matrix A stddev(A, norm type, dim) find the standard deviation in each column or row of A, using specified normalisation var(A, norm type, dim) find the variance in each column or row of matrix A, using specified normalisation   A.print("A:") print the entire contents of A save -ascii 'A.txt' A A.save("A.txt", raw ascii) Matlab/Octave matrices saved as ascii text are readable load -ascii 'A.txt' A.load("A.txt", raw ascii) by PyArmadillo (and vice-versa) A = rand(2,3) A = randu(2,3) randu generates uniformly distributed random numbers B = randn(4,5) B = randn (4,5)