jax.scipy.linalg.solve_sylvester#

jax.scipy.linalg.solve_sylvester(A, B, C, *, method='schur', tol=1e-08)[source]#

Solves the Sylvester equation .. math:

AX + XB = C

Using one of two methods.

  1. Bartell-Stewart (schur) algorithm (default) [CPU ONLY]:

Where A and B are first decomposed using Schur decomposition to construct and alternate sylvester equation: .. math:

RY + YS^T = F

Where R and S are in quasitriangular form when A and B are real valued and triangular when A and B are complex.

  1. The Eigen decomposition algorithm [CPU and GPU]

Parameters:

  • A (ArrayLike) – Matrix of shape m x m

  • B (ArrayLike) – Matrix of shape n x n

  • C (ArrayLike) – Matrix of shape m x n

  • method (str) – “schur” is the default and is accurate but slow, and “eigen” is an alternative that is faster but less accurate for ill-conditioned matrices.

  • tol (float) – How close the sum of the eigenvalues from A and B can be to zero before returning matrix of NaNs

Returns:

Matrix of shape m x n

Return type:

X

Examples

>>> A = jax.numpy.array([[1, 2], [3, 4]])
>>> B = jax.numpy.array([[5, 6], [7, 8]])
>>> C = jax.numpy.array([[6, 8], [10, 12]])
>>> X = jax.scipy.linalg.solve_sylvester(A, B, C)
>>> print(X) 
[[1. 0.]
 [0.  1.]]

Notes

The Bartel-Stewart algorithm is robust because a Schur decomposition always exists even for defective matrices, and it handles complex and ill-conditioned problems better than the eigen decomposition method. However, there are a couple of drawbacks. First, It is computationally more expensive than the eigen decomposition method because you need to perform a Schur decomposition and then scan the entire solution matrix. Second, it requires more system memory compared to the eigen decomposition method.

The eigen decomposition method is the fastest method to solve a sylvester equation. However, this speed brings with it a couple of drawbacks. First, A and B must be diagonalizable otherwise the eigenvectors will be linearly dependent and ill-conditioned leading to accuracy issues. Second, when the eigenvectors are not orthogonal roundoff errors are amplified.

Additionally, for complex types as the size of the matrix increases the accuracy of the results degrades. Float64 types are most robust to degradation.

The tol argument allows you to specify how ill-conditioned a matrix can be and still estimate a solution. For matrices that are ill-conditioned we recommend using float64 instead of the default float32 dtype. The solver can still return good estimates for ill-conditioned matrices depending on how close to zero the sums of the eigenvalues of A and B are.