Sparse iterative solvers (gmres, cg) in SciPy

Introduction

Sparse iterative solvers help solve big systems of equations quickly when most values are zero. They save time and memory.

When solving large linear systems from engineering or physics problems with many zero values.

When you want to find approximate solutions fast instead of exact ones.

When the matrix is too big to use direct methods like matrix inversion.

When working with sparse data in machine learning or network analysis.

When you want to save computer memory by avoiding dense matrix operations.

Syntax

SciPy

from scipy.sparse.linalg import gmres, cg

# gmres(A, b, tol=1e-5, maxiter=None)
# cg(A, b, tol=1e-5, maxiter=None)

# A: sparse matrix or linear operator
# b: right-hand side vector
# tol: tolerance for convergence
# maxiter: maximum iterations allowed

gmres works well for general matrices.

cg is faster but only works if the matrix is symmetric and positive definite.

Examples

Use gmres to solve Ax = b for a sparse matrix A.

SciPy

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

A = diags([2, 1, 3], [0, -1, 1], shape=(3,3))
b = np.array([1, 2, 3])
x, info = gmres(A, b)
print(x)

Use cg for symmetric positive definite sparse matrix A.

SciPy

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

A = diags([4, 1, 1], [0, -1, 1], shape=(3,3))
b = np.array([1, 2, 3])
x, info = cg(A, b)
print(x)

Sample Program

This program creates a small sparse matrix and solves the system Ax = b using both gmres and cg solvers. It prints the solutions and info codes where 0 means the solver succeeded.

SciPy

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

# Create a sparse matrix A
# Diagonal values: main diagonal 4, sub diagonal 1, super diagonal 1
A = diags([4, 1, 1], [0, -1, 1], shape=(3,3))

# Right-hand side vector b
b = np.array([1, 2, 3])

# Solve using gmres
x_gmres, info_gmres = gmres(A, b, tol=1e-8)

# Solve using cg
x_cg, info_cg = cg(A, b, tol=1e-8)

print("Solution with gmres:", x_gmres)
print("Info gmres (0 means success):", info_gmres)
print("Solution with cg:", x_cg)
print("Info cg (0 means success):", info_cg)

Important Notes

Check the info output: 0 means the solver found a solution successfully.

Set tol smaller for more accurate results but slower computation.

cg requires the matrix to be symmetric and positive definite; otherwise, it may fail.

Summary

Sparse iterative solvers like gmres and cg solve large sparse linear systems efficiently.

gmres works for general matrices; cg is faster but needs symmetric positive definite matrices.

They save memory and time by using the sparse structure instead of dense matrix methods.

Practice

1. Which of the following statements about scipy.sparse.linalg.cg is true?

easy

A. cg is slower than gmres for all matrices.

B. cg can solve any linear system regardless of matrix properties.

C. cg uses dense matrix methods internally.

D. cg requires the matrix to be symmetric and positive definite.

Solution

1. Step 1: Understand the requirements of cg

   The conjugate gradient method (cg) is designed for symmetric positive definite matrices only.

2. Step 2: Compare with other options

   cg cannot solve any matrix (B is wrong), it is often faster than gmres for suitable matrices (A is wrong), and it uses sparse methods, not dense (D is wrong).

3. Final Answer:

   cg requires the matrix to be symmetric and positive definite. -> Option D

4. Quick Check:

   cg needs symmetric positive definite matrix [OK]

Hint: Remember: CG needs symmetric positive definite matrices [OK]

Common Mistakes:

• Thinking CG works for any matrix
• Confusing CG with GMRES
• Assuming CG uses dense matrix methods

2. Which is the correct way to import the GMRES solver from SciPy?

easy

A. import scipy.linalg.gmres

B. from scipy.sparse.linalg import gmres

C. from scipy.linalg import gmres

D. import gmres from scipy.sparse

Solution

1. Step 1: Recall the module location of GMRES

   The GMRES solver is in scipy.sparse.linalg, so it must be imported from there.

2. Step 2: Check the import syntax

   Correct Python import syntax for a function is from module import function. from scipy.sparse.linalg import gmres matches this and the correct module.

3. Final Answer:

   from scipy.sparse.linalg import gmres -> Option B

4. Quick Check:

   Correct import syntax and module [OK]

Hint: Import gmres from scipy.sparse.linalg using 'from ... import' [OK]

Common Mistakes:

• Importing from scipy.linalg instead of sparse.linalg
• Using incorrect import syntax
• Trying to import gmres directly from scipy.sparse

3. What will be the output of the following code snippet?

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

A = diags([1, 2, 1], [-1, 0, 1]).toarray()
b = np.array([4, 6, 4])
x, info = cg(A, b)
print(np.round(x, 2))

medium

A. [1. 2. 1.]

B. [2. 1. 2.]

C. [0. 0. 0.]

D. Error due to matrix not positive definite

Solution

1. Step 1: Analyze the matrix and vector

   The matrix A is tridiagonal with diagonals [1,2,1], which is symmetric positive definite. Vector b is [4,6,4].

2. Step 2: Solve using conjugate gradient

   Using cg, the solution x satisfies Ax = b. The solution is approximately [1, 2, 1].

3. Final Answer:

   [1. 2. 1.] -> Option A

4. Quick Check:

   cg solves Ax=b with symmetric positive definite A [OK]

Hint: Check matrix symmetry and positive definiteness before cg [OK]

Common Mistakes:

• Assuming cg fails on this matrix
• Confusing the solution vector with b
• Not rounding output before comparing

4. Identify the error in this code using cg solver:

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

A = np.array([[0, 1], [1, 0]])
b = np.array([1, 2])
x, info = cg(A, b)
print(x)

medium

A. Matrix A is not symmetric positive definite, so cg will fail.

B. Vector b has wrong shape.

C. cg function is not imported correctly.

D. No error; code runs fine.

Solution

1. Step 1: Check matrix properties

   Matrix A = [[0,1],[1,0]] is symmetric but not positive definite (its eigenvalues are 1 and -1).

2. Step 2: Understand cg requirements

   The conjugate gradient method requires A to be symmetric positive definite. Since A is not, cg will fail or not converge properly.

3. Final Answer:

   Matrix A is not symmetric positive definite, so cg will fail. -> Option A

4. Quick Check:

   cg needs symmetric positive definite matrix [OK]

Hint: Check matrix eigenvalues before using cg [OK]

Common Mistakes:

• Ignoring matrix definiteness
• Assuming cg works for any symmetric matrix
• Thinking vector shape causes error

5. You have a large sparse matrix that is not symmetric positive definite. Which solver should you use to solve Ax = b efficiently?

hard

A. Use cg because it is always faster.

B. Convert matrix to dense and use direct solver.

C. Use gmres because it works for general matrices.

D. Use cg after transposing the matrix.

Solution

1. Step 1: Identify matrix properties

   The matrix is large, sparse, and not symmetric positive definite, so cg is not suitable.

2. Step 2: Choose appropriate solver

   gmres can handle general matrices efficiently without requiring symmetry or positive definiteness.

3. Step 3: Avoid dense conversion

   Converting to dense wastes memory and time, so it is not efficient.

4. Final Answer:

   Use gmres because it works for general matrices. -> Option C

5. Quick Check:

   gmres handles general sparse matrices [OK]

Hint: Use gmres for non-symmetric or indefinite sparse matrices [OK]

Common Mistakes:

• Trying to use cg on non-symmetric matrices
• Converting sparse to dense unnecessarily
• Thinking transposing fixes definiteness