I am using scipy.sparse.linalg.spsolve() to solve a set of linear equations in a loop where the equations are opdated iteratively based on the previous solution. Sometimes I end up in a singular matrix problem where I use a try-except to add terms on the diagonal (Tikhonov regularization) before trying to solve again. This seem to work for the warning about singular matrix, but now I have reached another problem where it crashes in a way I do not catch with the try-except. The message is the following where I replaced actual names and paths on my system with "...":
Windows fatal exception: access violation
Main thread:
Current thread 0x00000270 (most recent call first):
File "...\lib\site-packages\scipy\sparse\linalg\dsolve\linsolve.py", line 203 in spsolve
File "...\curvefit.py", line 401 in solve_linear_system_with_regularization
File "...\curvefit.py", line 320 in smoothing_spline
File "...\curvefit.py", line 580 in run_curvefit
File "...\fitSplinesShapeConstraints.py", line 53 in fitSplinesShapeConstraints
File "...\main.py", line 107 in <module>
File "C:\Users\...\AppData\Local\anaconda3\envs\...\lib\site-packages\spyder_kernels\py3compat.py", line 356 in compat_exec
File "C:\Users\...\AppData\Local\anaconda3\envs\...\lib\site-packages\spyder_kernels\customize\spydercustomize.py", line 473 in exec_code
File "C:\Users\...\AppData\Local\anaconda3\envs\...\lib\site-packages\spyder_kernels\customize\spydercustomize.py", line 615 in _exec_file
File "C:\Users\...\AppData\Local\anaconda3\envs\...\lib\site-packages\spyder_kernels\customize\spydercustomize.py", line 528 in runfile
File "C:\Users\...\AppData\Local\Temp\ipykernel_25612\3794581249.py", line 1 in <module>
Restarting kernel...
My code looks like:
import numpy as np
import scipy.sparse as sp
import scipy.sparse.linalg as la
#from scipy.linalg import LinAlgError
from scipy.optimize import minimize_scalar
from scipy.sparse import identity
import scipy.sparse
from scipy.sparse.linalg.dsolve.linsolve import MatrixRankWarning # Import the custom warning class
import warnings
def solve_linear_system_with_regularization(A, bb):
regularization_range = np.logspace(-20, 15, 36) # Regularization parameter range
best_residual_norm = np.inf
best_solution = None
best_reg_param = None
consecutive_error_increases_regularized_solution = 0
with warnings.catch_warnings():
# Convert warnings to exceptions
warnings.simplefilter("error", MatrixRankWarning)
try:
parameters = la.spsolve(A, bb)
return parameters
except MatrixRankWarning as e:
# If the matrix is singular, apply Tikhonov regularization
if 'singular' in str(e):
print("Matrix is singular. Applying Tikhonov regularization in an attempt to solve anyways...")
for reg_param in regularization_range:
try:
print('trying with regularization parameter: ' + str(reg_param))
# Try solving the regularized system
# Create a sparse identity matrix with the same shape as A
identity_matrix = identity(A.shape[0], format='csr')
A_reg = A + reg_param * identity_matrix
x_reg = la.spsolve(A_reg, bb) # Regularized solution
# Calculate the dot product result
dot_product_result = np.dot(A.toarray(), x_reg)
# Calculate the residual norm (error) of the solution
residual_norm = np.linalg.norm(dot_product_result - bb)
# Check if the residual norm is increasing
if residual_norm > best_residual_norm:
consecutive_error_increases_regularized_solution += 1
else:
consecutive_error_increases_regularized_solution = 0
# Update best solution if the residual norm is smaller
if residual_norm < best_residual_norm:
best_solution = x_reg
best_reg_param = reg_param
best_residual_norm = residual_norm
except MatrixRankWarning:
# Skip if regularization causes singularity
continue
# Break out of the loop if residual_norm is increasing three times
if consecutive_error_increases_regularized_solution >= 3:
print('The error from adding more regularization seems to be increasing for the last 3 increases of the regularization parameter')
break
if best_solution is None:
raise Warning("Unable to find a suitable regularization parameter.")
print("The system of equations has been solved with a regularization parameter of:", best_reg_param)
return best_solution
else:
# If it's another type of error, raise it
# I never reach this point...
raise e
I have thought about adding other things for the except line to catch like except (MatrixRankWarning, ??) as e:
but I do not know what to put at the '??' to catch this "Windows fatal exception: access violation"
A solution to avoiding the error (not catching it as this is not possible) is:
scipy.sparse.linalg..spsolve(A_reg, bb, permc_spec='NATURAL')
instead of the default:
permc_spec='COLAMD'.
I do not know why they perform differently in this case.