def regular_inv(A, EPS_A=1e-07, EPS_residual=0.001):
    """Computes the inverse of a matrix. Attempts a regularization by adding an epsilon to the diagonal, if the matrix is ill-conditioned.

    Parameters
    ----------
    A : array-like, square
        input matrix
    EPS_A : double, optional, default=1e-7
        epsilon to be added to the diagonal of A, if necessary
    EPS_residual : double, optional, default=1e-3
        maximum residual in A*A^{-1} before attempting regularization

    Returns
    -------
    A^{-1} : array-like, square
        numerical inverse of A

    """
    import numpy as np
    import scipy.linalg as sla
    inv_A = sla.inv(A)
    residual = np.sum(np.fabs(inv_A.dot(A) - np.identity(A.shape[0])))
    if residual > EPS_residual:
        A = A + EPS_A * np.identity(A.shape[0])
        inv_A = sla.inv(A)
        residual = np.sum(np.fabs(inv_A.dot(A) - np.identity(A.shape[0])))
        if residual > 0.001:
            from pyselfi.utils import PrintValue
            PrintValue('condition number of A', np.linalg.cond(A))
            PrintValue('residuals in A*A^{-1}', residual)
            raise ValueError(
                'Matrix inversion is imprecise for A, even after regularization attempt.'
                )
        inv_A = inv_A - EPS_A * np.identity(A.shape[0])
    return inv_A
