linalg._matfuncs_inv_ssq

Matrix functions that use Pade approximation with inverse scaling and squaring.

Module Contents

Classes

LogmRankWarning()
LogmExactlySingularWarning()
LogmNearlySingularWarning()
LogmError()
FractionalMatrixPowerError()
_MatrixM1PowerOperator(self,A,p) A representation of the linear operator (A - I)^p.

Functions

_onenormest_m1_power(A,p,t=2,itmax=5,compute_v=False,compute_w=False) Efficiently estimate the 1-norm of (A - I)^p.
_unwindk(z) Compute the scalar unwinding number.
_briggs_helper_function(a,k) Computes r = a^(1 / (2^k)) - 1.
_fractional_power_superdiag_entry(l1,l2,t12,p) Compute a superdiagonal entry of a fractional matrix power.
_logm_superdiag_entry(l1,l2,t12) Compute a superdiagonal entry of a matrix logarithm.
_inverse_squaring_helper(T0,theta) A helper function for inverse scaling and squaring for Pade approximation.
_fractional_power_pade_constant(i,t)
_fractional_power_pade(R,t,m) Evaluate the Pade approximation of a fractional matrix power.
_remainder_matrix_power_triu(T,t) Compute a fractional power of an upper triangular matrix.
_remainder_matrix_power(A,t) Compute the fractional power of a matrix, for fractions -1 < t < 1.
_fractional_matrix_power(A,p) Compute the fractional power of a matrix.
_logm_triu(T) Compute matrix logarithm of an upper triangular matrix.
_logm_force_nonsingular_triangular_matrix(T,inplace=False)
_logm(A) Compute the matrix logarithm.
class LogmRankWarning
class LogmExactlySingularWarning
class LogmNearlySingularWarning
class LogmError
class FractionalMatrixPowerError
class _MatrixM1PowerOperator(A, p)

A representation of the linear operator (A - I)^p.

__init__(A, p)
_matvec(x)
_rmatvec(x)
_matmat(X)
_adjoint()
_onenormest_m1_power(A, p, t=2, itmax=5, compute_v=False, compute_w=False)

Efficiently estimate the 1-norm of (A - I)^p.

A : ndarray
Matrix whose 1-norm of a power is to be computed.
p : int
Non-negative integer power.
t : int, optional
A positive parameter controlling the tradeoff between accuracy versus time and memory usage. Larger values take longer and use more memory but give more accurate output.
itmax : int, optional
Use at most this many iterations.
compute_v : bool, optional
Request a norm-maximizing linear operator input vector if True.
compute_w : bool, optional
Request a norm-maximizing linear operator output vector if True.
est : float
An underestimate of the 1-norm of the sparse matrix.
v : ndarray, optional
The vector such that ||Av||_1 == est*||v||_1. It can be thought of as an input to the linear operator that gives an output with particularly large norm.
w : ndarray, optional
The vector Av which has relatively large 1-norm. It can be thought of as an output of the linear operator that is relatively large in norm compared to the input.
_unwindk(z)

Compute the scalar unwinding number.

Uses Eq. (5.3) in [1]_, and should be equal to (z - log(exp(z)) / (2 pi i). Note that this definition differs in sign from the original definition in equations (5, 6) in [2]_. The sign convention is justified in [3]_.

z : complex
A complex number.
unwinding_number : integer
The scalar unwinding number of z.
[1]Nicholas J. Higham and Lijing lin (2011) “A Schur-Pade Algorithm for Fractional Powers of a Matrix.” SIAM Journal on Matrix Analysis and Applications, 32 (3). pp. 1056-1078. ISSN 0895-4798
[2]Robert M. Corless and David J. Jeffrey, “The unwinding number.” Newsletter ACM SIGSAM Bulletin Volume 30, Issue 2, June 1996, Pages 28-35.
[3]Russell Bradford and Robert M. Corless and James H. Davenport and David J. Jeffrey and Stephen M. Watt, “Reasoning about the elementary functions of complex analysis” Annals of Mathematics and Artificial Intelligence, 36: 303-318, 2002.
_briggs_helper_function(a, k)

Computes r = a^(1 / (2^k)) - 1.

This is algorithm (2) of [1]_. The purpose is to avoid a danger of subtractive cancellation. For more computational efficiency it should probably be cythonized.

a : complex
A complex number.
k : integer
A nonnegative integer.
r : complex
The value r = a^(1 / (2^k)) - 1 computed with less cancellation.

The algorithm as formulated in the reference does not handle k=0 or k=1 correctly, so these are special-cased in this implementation. This function is intended to not allow a to belong to the closed negative real axis, but this constraint is relaxed.

[1]Awad H. Al-Mohy (2012) “A more accurate Briggs method for the logarithm”, Numerical Algorithms, 59 : 393–402.
_fractional_power_superdiag_entry(l1, l2, t12, p)

Compute a superdiagonal entry of a fractional matrix power.

This is Eq. (5.6) in [1]_.

l1 : complex
A diagonal entry of the matrix.
l2 : complex
A diagonal entry of the matrix.
t12 : complex
A superdiagonal entry of the matrix.
p : float
A fractional power.
f12 : complex
A superdiagonal entry of the fractional matrix power.

Care has been taken to return a real number if possible when all of the inputs are real numbers.

[1]Nicholas J. Higham and Lijing lin (2011) “A Schur-Pade Algorithm for Fractional Powers of a Matrix.” SIAM Journal on Matrix Analysis and Applications, 32 (3). pp. 1056-1078. ISSN 0895-4798
_logm_superdiag_entry(l1, l2, t12)

Compute a superdiagonal entry of a matrix logarithm.

This is like Eq. (11.28) in [1]_, except the determination of whether l1 and l2 are sufficiently far apart has been modified.

l1 : complex
A diagonal entry of the matrix.
l2 : complex
A diagonal entry of the matrix.
t12 : complex
A superdiagonal entry of the matrix.
f12 : complex
A superdiagonal entry of the matrix logarithm.

Care has been taken to return a real number if possible when all of the inputs are real numbers.

[1]Nicholas J. Higham (2008) “Functions of Matrices: Theory and Computation” ISBN 978-0-898716-46-7
_inverse_squaring_helper(T0, theta)

A helper function for inverse scaling and squaring for Pade approximation.

T0 : (N, N) array_like upper triangular
Matrix involved in inverse scaling and squaring.
theta : indexable
The values theta[1] .. theta[7] must be available. They represent bounds related to Pade approximation, and they depend on the matrix function which is being computed. For example, different values of theta are required for matrix logarithm than for fractional matrix power.
R : (N, N) array_like upper triangular
Composition of zero or more matrix square roots of T0, minus I.
s : non-negative integer
Number of square roots taken.
m : positive integer
The degree of the Pade approximation.

This subroutine appears as a chunk of lines within a couple of published algorithms; for example it appears as lines 4–35 in algorithm (3.1) of [1]_, and as lines 3–34 in algorithm (4.1) of [2]_. The instances of ‘goto line 38’ in algorithm (3.1) of [1]_ probably mean ‘goto line 36’ and have been intepreted accordingly.

[1]Nicholas J. Higham and Lijing Lin (2013) “An Improved Schur-Pade Algorithm for Fractional Powers of a Matrix and their Frechet Derivatives.”
[2]Awad H. Al-Mohy and Nicholas J. Higham (2012) “Improved Inverse Scaling and Squaring Algorithms for the Matrix Logarithm.” SIAM Journal on Scientific Computing, 34 (4). C152-C169. ISSN 1095-7197
_fractional_power_pade_constant(i, t)
_fractional_power_pade(R, t, m)

Evaluate the Pade approximation of a fractional matrix power.

Evaluate the degree-m Pade approximation of R to the fractional matrix power t using the continued fraction in bottom-up fashion using algorithm (4.1) in [1]_.

R : (N, N) array_like
Upper triangular matrix whose fractional power to evaluate.
t : float
Fractional power between -1 and 1 exclusive.
m : positive integer
Degree of Pade approximation.
U : (N, N) array_like
The degree-m Pade approximation of R to the fractional power t. This matrix will be upper triangular.
[1]Nicholas J. Higham and Lijing lin (2011) “A Schur-Pade Algorithm for Fractional Powers of a Matrix.” SIAM Journal on Matrix Analysis and Applications, 32 (3). pp. 1056-1078. ISSN 0895-4798
_remainder_matrix_power_triu(T, t)

Compute a fractional power of an upper triangular matrix.

The fractional power is restricted to fractions -1 < t < 1. This uses algorithm (3.1) of [1]_. The Pade approximation itself uses algorithm (4.1) of [2]_.

T : (N, N) array_like
Upper triangular matrix whose fractional power to evaluate.
t : float
Fractional power between -1 and 1 exclusive.
X : (N, N) array_like
The fractional power of the matrix.
[1]Nicholas J. Higham and Lijing Lin (2013) “An Improved Schur-Pade Algorithm for Fractional Powers of a Matrix and their Frechet Derivatives.”
[2]Nicholas J. Higham and Lijing lin (2011) “A Schur-Pade Algorithm for Fractional Powers of a Matrix.” SIAM Journal on Matrix Analysis and Applications, 32 (3). pp. 1056-1078. ISSN 0895-4798
_remainder_matrix_power(A, t)

Compute the fractional power of a matrix, for fractions -1 < t < 1.

This uses algorithm (3.1) of [1]_. The Pade approximation itself uses algorithm (4.1) of [2]_.

A : (N, N) array_like
Matrix whose fractional power to evaluate.
t : float
Fractional power between -1 and 1 exclusive.
X : (N, N) array_like
The fractional power of the matrix.
[1]Nicholas J. Higham and Lijing Lin (2013) “An Improved Schur-Pade Algorithm for Fractional Powers of a Matrix and their Frechet Derivatives.”
[2]Nicholas J. Higham and Lijing lin (2011) “A Schur-Pade Algorithm for Fractional Powers of a Matrix.” SIAM Journal on Matrix Analysis and Applications, 32 (3). pp. 1056-1078. ISSN 0895-4798
_fractional_matrix_power(A, p)

Compute the fractional power of a matrix.

See the fractional_matrix_power docstring in matfuncs.py for more info.

_logm_triu(T)

Compute matrix logarithm of an upper triangular matrix.

The matrix logarithm is the inverse of expm: expm(logm(T)) == T

T : (N, N) array_like
Upper triangular matrix whose logarithm to evaluate
logm : (N, N) ndarray
Matrix logarithm of T
[1]Awad H. Al-Mohy and Nicholas J. Higham (2012) “Improved Inverse Scaling and Squaring Algorithms for the Matrix Logarithm.” SIAM Journal on Scientific Computing, 34 (4). C152-C169. ISSN 1095-7197
[2]Nicholas J. Higham (2008) “Functions of Matrices: Theory and Computation” ISBN 978-0-898716-46-7
[3]Nicholas J. Higham and Lijing lin (2011) “A Schur-Pade Algorithm for Fractional Powers of a Matrix.” SIAM Journal on Matrix Analysis and Applications, 32 (3). pp. 1056-1078. ISSN 0895-4798
_logm_force_nonsingular_triangular_matrix(T, inplace=False)
_logm(A)

Compute the matrix logarithm.

See the logm docstring in matfuncs.py for more info.

In this function we look at triangular matrices that are similar to the input matrix. If any diagonal entry of such a triangular matrix is exactly zero then the original matrix is singular. The matrix logarithm does not exist for such matrices, but in such cases we will pretend that the diagonal entries that are zero are actually slightly positive by an ad-hoc amount, in the interest of returning something more useful than NaN. This will cause a warning.