pythonscipytry-except

How can I catch exeption of scipy.sparse.linalg.spsolve() (Windows fatal exception: access violation)


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"


Solution

  • 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.