Preconditioners in SciPy - Mini Project: Build & Apply

Using Preconditioners to Solve Linear Systems Efficiently

📖 Scenario: You work as a data scientist helping engineers solve large systems of linear equations quickly. These systems come from real-world problems like simulating airflow or electrical circuits.Sometimes, solving these systems directly takes too long. Preconditioners are tools that help speed up the solving process by improving the system's properties.

🎯 Goal: Build a Python program that creates a linear system, sets up a simple preconditioner, solves the system using an iterative solver with the preconditioner, and prints the solution.

📋 What You'll Learn

Create a sparse matrix representing the system

Create a right-hand side vector

Create a Jacobi preconditioner

Use the conjugate gradient solver with the preconditioner

Print the solution vector

💡 Why This Matters

🌍 Real World

Preconditioners are used in engineering and scientific computing to speed up solving large systems of equations from simulations and models.

💼 Career

Understanding preconditioners helps data scientists and engineers optimize numerical algorithms for faster and more stable solutions.

1

Create the linear system matrix and vector

Create a sparse matrix called A using scipy.sparse.diags with diagonals [4, -1, -1] and offsets [0, -1, 1] for size 5. Create a vector b as a NumPy array with values [1, 2, 3, 4, 5].

SciPy

# Create the sparse matrix A and vector b
# Your code here

Need a hint?

Use diags to create a tridiagonal matrix with main diagonal 4 and -1 on the diagonals just above and below.

2

Create the Jacobi preconditioner

Create a diagonal preconditioner matrix M_inv as the inverse of the diagonal of A using scipy.sparse.diags. Extract the diagonal of A with A.diagonal() and invert it element-wise.

SciPy

import numpy as np
from scipy.sparse import diags

A = diags([4, -1, -1], [0, -1, 1], shape=(5, 5))
b = np.array([1, 2, 3, 4, 5])

# Create the Jacobi preconditioner M_inv
# Your code here

Need a hint?

Get the diagonal of A with A.diagonal(). Then create M_inv with the inverse values on the diagonal.

3

Solve the system using conjugate gradient with preconditioner

Use scipy.sparse.linalg.cg to solve Ax = b with the preconditioner M_inv. Pass M_inv as the M argument. Store the solution vector in x.

SciPy

import numpy as np
from scipy.sparse import diags
from scipy.sparse.linalg import cg

A = diags([4, -1, -1], [0, -1, 1], shape=(5, 5))
b = np.array([1, 2, 3, 4, 5])

diag_A = A.diagonal()
M_inv = diags(1 / diag_A)

# Solve Ax = b using conjugate gradient with preconditioner M_inv
# Your code here

Need a hint?

Call cg with A, b, and M=M_inv. Capture the solution in x.

4

Print the solution vector

Print the solution vector x using print(x).

SciPy

import numpy as np
from scipy.sparse import diags
from scipy.sparse.linalg import cg

A = diags([4, -1, -1], [0, -1, 1], shape=(5, 5))
b = np.array([1, 2, 3, 4, 5])

diag_A = A.diagonal()
M_inv = diags(1 / diag_A)

x, info = cg(A, b, M=M_inv)

# Print the solution vector x
# Your code here

Need a hint?

Use print(x) to display the solution vector.