SciPy - solve_sylvester() Function

---

---

The solve_sylvester() function in SciPy aims to find a solution to the Sylvester equation AX+XB=Q. In this equation, A B, and Q are given matrices, while X stands for the solution matrix we need to find.

This function relies on a number-based approach that builds on the Bartels-Stewart algorithm. This algorithm changes matrices A and B into triangle shapes through Schur decomposition. This change makes the Sylvester equation easier to solve. The transformation ensures stability and good performance even with large matrices. Many experts use this method in computer-based linear algebra.

Engineers often use the solve_sylvester method to check if a system will reach a steady state. This method solves the Lyapunov equation, which is a type of Sylvester's equation. The solution gives us a matrix X. The system's stability depends on this matrix's eigen values. Positive eigen values mean the system is stable and will settle down over time. Non-positive definite eigen values indicate an unstable or marginally stable.

Knowing about stability is important because it helps ensure the systems like autopilots or factory control setups work and as expected in the long run.This equation plays a key role in control theory, signal processing, and other advanced applications of linear algebra.

Syntax

The syntax for the Scipy solve_sylvester() method is as follows −

.solve_sylvester(a, b, q)

Parameters

This method accepts the following parameters −

• a array like shape (m,m) − Coefficient matrix A

• b array like shape (n,n)− Coefficient matrix B

• q array like shape (m,n)− Matrix Q on the right-hand side of the Sylvester equation.

Return Value

x ndarray(m,n) − The solution matrix X that satisfies AX+XB=Q.

Example 1: Solving Sylvester Equation

This is the basic example of SciPy solve_sylvester() method where we have created three matrices A, B, Q to solve the X for the Sylvester equation AX=XB=Q.

import numpy as np
from scipy import linalg
# Define matrices A, B, and Q
A = np.array([[1, 2], [3, 4]])
B = np.array([[2, 0], [0, -3]])
Q = np.array([[1, 0], [0, 1]])

# Solve the Sylvester equation
X = linalg.solve_sylvester(A, B, Q)

print("Solution Matrix X:\n", X)

When we run above program, it produces following result

Solution Matrix X:
 [[ 0.5   0.25]
 [-0.25  0.25]]

Example 2: Stability Analysis Using Lyapunov Equation

In the below example we will solve the specific form of Sylvester equation i.e, Lyapunov equation AX+XA^-T=-Q using solve_sylvester() method. This is used to determine the stability of the system.

import numpy as np
from scipy.linalg import solve_sylvester, eigvals

A = np.array([[0, 1], [-2, -3]])
Q = np.eye(2)  # Identity matrix

# Solve the Lyapunov equation AX + XA^T = -Q
X = solve_sylvester(A, A.T, -Q)

# Check if X is positive definite by verifying eigenvalues
eigenvalues = eigvals(X)

print("Solution to Lyapunov Equation (X):\n", X)
print("\nEigenvalues of X:", eigenvalues)
print("\nIs X Positive Definite?:", np.all(eigenvalues > 0))

Following is an output of the above code

Solution to Lyapunov Equation (X):
 [[ 1.  -0.5]
 [-0.5  0.5]]

Eigenvalues of X: [1.30901699+0.j 0.19098301+0.j]

Is X Positive Definite?: True

Example 3: Combining Lyapunov Equation and expm() method

In the below code we will solve the Lyapunov equation using solve_sylvester() method. We the use the expm() method to compute the system dynamic over time, the matrix 'A' helps the model to know how the state 'X' of the system evolves over time and we use solution matrix 'X' to find the stability of the system. Following is the code −

import numpy as np
from scipy.linalg import solve_sylvester, expm, eig

A = np.array([[0, 1], [-2, -3]])
Q = np.eye(2)  # Identity matrix
X = solve_sylvester(A, A.T, -Q)
print("Solution to Lyapunov Equation (X):\n", X)

# Combine with matrix exponential to study system dynamics over time
time = 5
system_dynamics = expm(A * time)
print("\nSystem dynamics at time t =", time, ":\n", system_dynamics)

# Eigenvalue analysis for stability conditions
eigenvalues = eigvals(X)
print("\nEigenvalues of X:", eigenvalues)
print("\nIs X Positive Definite?:", np.all(eigenvalues > 0))

Output of the above code is as follows

Solution to Lyapunov Equation (X):
 [[ 1.  -0.5]
 [-0.5  0.5]]

System dynamics at time t = 5 :
 [[ 0.01343049  0.00669255]
 [-0.01338509 -0.00664715]]

Eigenvalues of X: [1.30901699+0.j 0.19098301+0.j]

Is X Positive Definite?: True