A matrix function is a function that maps a square matrix to another square matrix via a power series. These should not be confused with vectorization: the application of any given function of one variable to each element of a matrix. For example, it is not the same to compute the square of a square matrix, A.dot(A)
(for example, In [8]
), than a matrix with all the elements of A
squared (examples In [5]
through In []
).
Note
To make the proper distinction in notation, we will write A^2
to denote the actual square of a square matrix and A^n
to represent the subsequent powers (for all positive integers n).
Constructor 
Description 


Extension of a scalarvalued function called 

Fractional matrix power 

Matrix exponential 

Action of the matrix exponential of 

Frechet derivative of the matrix exponential in the 

Matrix cosine 

Matrix sine 

Matrix tangent 

Hyperbolic matrix cosine 

Hyperbolic matrix sine 

Hyperbolic matrix tangent 

Matrix sign function 

Matrix square root 

Matrix logarithm 
In [1]: import numpy as np, scipy as sp; \ ...: import scipy.linalg as spla In [2]: np.set_printoptions(suppress=True, precision=3) In [3]: def square(x): return x**2 In [4]: A = spla.hilbert(4); print A [[ 1. 0.5 0.333 0.25 ] [ 0.5 0.333 0.25 0.2 ] [ 0.333 0.25 0.2 0.167] [ 0.25 0.2 0.167 0.143]] In [5]: print square(A) [[ 1. 0.25 0.111 0.062] [ 0.5 0.333 0.25 0.2 ] [ 0.333 0.25 0.2 0.167] [ 0.25 0.2 0.167 0.143]] In [6]: print A*A [[ 1. 0.25 0.111 0.062] [ 0.25 0.111 0.062 0.04 ] [ 0.111 0.062 0.04 0.028] [ 0.062 0.04 0.028 0.02 ]] In [7]: print A**2 [[ 1. 0.25 0.111 0.062] [ 0.25 0.111 0.062 0.04 ] [ 0.111 0.062 0.04 0.028] [ 0.062 0.04 0.028 0.02 ]] In [8]: print A.dot(A) [[ 1.424 0.8 0.567 0.441] [ 0.8 0.464 0.333 0.262] [ 0.567 0.333 0.241 0.19 ] [ 0.441 0.262 0.19 0.151]]
The actual powers A^n
of a matrix is the starting point for the definition of any matrix function. In the module numpy.linalg
we have the routine matrix_power
to perform this operation. We can also achieve this result with the generic function funm
or with the function fractional_matrix_power
, both of them in the module scipy.linalg
.
In [9]: print np.linalg.matrix_power(A, 2) [[ 1.424 0.8 0.567 0.441] [ 0.8 0.464 0.333 0.262] [ 0.567 0.333 0.241 0.19 ] [ 0.441 0.262 0.19 0.151]] In [10]: print spla.fractional_matrix_power(A, 2) [[ 1.424 0.8 0.567 0.441] [ 0.8 0.464 0.333 0.262] [ 0.567 0.333 0.241 0.19 ] [ 0.441 0.262 0.19 0.151]] In [11]: print spla.funm(A, square) [[ 1.424 0.8 0.567 0.441] [ 0.8 0.464 0.333 0.262] [ 0.567 0.333 0.241 0.19 ] [ 0.441 0.262 0.19 0.151]]
To compute any matrix function, theoretically, we first express the function as a power series, by means of its Taylor expansion. Then, we apply the input matrix into an approximation to that expansion (since it is impossible to add matrices ad infinitum). Most matrix functions necessarily carry an error of computation, for this reason. In the scipy.linalg
module, the matrix functions are coded following this principle.
Note that there are three functions with an optional Boolean parameter
disp
. To understand the usage of this parameter, we must remember that most matrix functions compute approximations, with an error of computation. The parameterdisp
is set toTrue
by default, and it produces a warning if the error of approximation is large. If we setdisp
toFalse
, instead of a warning we will obtain the 1norm of the estimated error.The algorithms behind the functions
expm
, the action of an exponential over a matrix,expm_multiply
, and the Frechet derivative of an exponential,expm_frechet
, use Pade approximations instead of Taylor expansions. This allows for more robust and accurate calculations. All the trigonometric and hyperbolic trigonometric functions base their algorithm in easy computations involvingexpm
.The generic matrix function called
funm
and the squareroot function calledsqrtm
apply clever algorithms that play with the Schur decomposition of the input matrix, and proper algebraic manipulations with the corresponding eigenvalues. They are still prone to roundoff errors but are much faster and more accurate than any algorithm based on Taylor expansions.The matrix sign function called
signm
is initially an application offunm
with the appropriate function, but should this approach fail, the algorithm takes a different approach based on iterations that converges to a decent approximation to the solution.The functions
logm
andfractional_matrix_power
(when the latter is applied to noninteger powers) use a very complex combination (and improvement!) of Pade approximations and Schur decompositions.
Tip
We will explore Schur decompositions when we deal with matrix factorizations related to eigenvalues. In the meantime, if you are interested in learning the particulars of these clever algorithms, read their descriptions in Golub and Van Loan, Matrix Computations 4 edition, Johns Hopkins Studies in the Mathematical Sciences, vol. 3.
For details on the improvements to SchurPade algorithms, as well as the algorithm behind Frechet derivatives of the exponential, refer to:
Nicholas J. Higham and Lijing Lin An Improved SchurPade Algorithm for Fractional Powers of a Matrix and Their Frechet Derivatives
Awad H. AlMohy and Nicholas J. Higham Improved Inverse Scaling and Squaring Algorithms for the Matrix Logarithm, in SIAM Journal on Scientific Computing, 34 (4)