Sparse iterative solvers (gmres, cg) in SciPy - Practice Problems & Coding Challenges

Challenge - 5 Problems

❓ Predict Output

intermediate

Output of CG solver on a simple sparse system

What is the output of the following code snippet that uses the Conjugate Gradient (CG) solver from scipy.sparse.linalg to solve a sparse linear system?

SciPy

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

# Create a sparse diagonal matrix
n = 5
k = [-1*np.ones(n-1), 2*np.ones(n), -1*np.ones(n-1)]
offsets = [-1, 0, 1]
A = diags(k, offsets)

b = np.array([1, 2, 3, 4, 5])

x, info = cg(A, b)
print(x)

A [5.83 10.67 13.5 13.33 9.17]

B [2. 3. 4. 5. 6.]

C [0.75 1.5 2. 2.5 3. ]

D [1. 2. 3. 4. 5.]

💡 Hint

Think about how the CG solver solves Ax = b for a tridiagonal matrix with 2 on the diagonal and -1 on the off-diagonals.

The matrix A is a tridiagonal matrix representing a simple finite difference operator. The CG solver finds the solution x such that Ax = b. The output vector x is approximately [0.75, 1.5, 2.0, 2.5, 3.0].

❓ Predict Output

intermediate

GMRES solver output for a sparse system

What is the output of the following code that uses GMRES to solve a sparse linear system?

SciPy

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

n = 4
k = [np.ones(n), 2*np.ones(n), np.ones(n)]
offsets = [-1, 0, 1]
A = diags(k, offsets)

b = np.array([1, 2, 3, 4])

x, info = gmres(A, b)
print(x)

A [0.25 0.5 0.75 1. ]

B [0. 1. 0. 2. ]

C [1. 2. 3. 4.]

D [0.5 1. 1.5 2. ]

💡 Hint

GMRES solves Ax = b for a matrix with 2 on the diagonal and 1 on the off-diagonals.

The matrix A has 2 on the diagonal and 1 on the off-diagonals. The GMRES solver finds x such that Ax = b. The solution vector x is approximately [0.25, 0.5, 0.75, 1.0].

❓ data_output

advanced

Number of iterations for CG solver convergence

Given the following code that solves a sparse system using CG with a tolerance of 1e-8, what is the number of iterations taken for convergence?

SciPy

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

n = 10
k = [-1*np.ones(n-1), 2*np.ones(n), -1*np.ones(n-1)]
offsets = [-1, 0, 1]
A = diags(k, offsets)

b = np.arange(1, n+1)

x, info = cg(A, b, tol=1e-8, maxiter=1000)
print(info)

A 0

B 10

C 1000

D 5

💡 Hint

The info output from cg returns 0 if convergence is reached successfully.

The CG solver converges successfully and info returns 0 indicating successful convergence within the allowed iterations.

🔧 Debug

advanced

Identify the error in GMRES solver usage

What error will the following code raise when trying to solve a sparse system with GMRES?

SciPy

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

n = 3
k = [np.ones(n), 2*np.ones(n), np.ones(n)]
offsets = [-1, 0, 1]
A = diags(k, offsets)

b = np.array([1, 2])

x, info = gmres(A, b)
print(x)

A IndexError: index out of range

B ValueError: dimension mismatch between A and b

C TypeError: unsupported operand type(s) for +: 'int' and 'str'

D No error, prints solution vector

💡 Hint

Check the size of matrix A and vector b.

The matrix A is 3x3 but vector b has length 2, causing a dimension mismatch error when calling gmres.

🚀 Application

expert

Choosing solver based on matrix properties

You have a large sparse symmetric positive definite matrix A and vector b. Which solver is the best choice to efficiently solve Ax = b?

A Use GMRES solver because it works for any matrix

B Use Jacobi iterative method for faster convergence

C Use direct solver like LU decomposition for sparse matrices

D Use CG solver because it is optimized for symmetric positive definite matrices

💡 Hint

Consider matrix properties and solver strengths.

CG solver is designed for symmetric positive definite matrices and is more efficient than GMRES or Jacobi for such cases. Direct solvers are often slower for large sparse matrices.