ill conditioned matrix in solving Ax=b

import numpy as np

from scipy.stats import ortho_group
from scipy.sparse import diags
from scipy import linalg

m = ortho_group.rvs(dim=5)
n = ortho_group.rvs(dim=5)

diags([5,4,3,2,1], 0).toarray()

array([[5., 0., 0., 0., 0.],
       [0., 4., 0., 0., 0.],
       [0., 0., 3., 0., 0.],
       [0., 0., 0., 2., 0.],
       [0., 0., 0., 0., 1.]])

A_2(A_4?) is ill-conditioned matrix

A_1 =  m.dot(diags([5,4,3,2,1], 0).toarray()).dot(n)
A_2 =  m.dot(diags([5,4,3,2,0.0001], 0).toarray()).dot(n)
A_3 =  m.dot(diags([5,4,3,2,0], 0).toarray()).dot(n)
A_4 =  m.dot(diags([5,4,3,2,0.1], 0).toarray()).dot(n)

A_1, A_2, A_3, A_4

(array([[-0.55566913,  2.04645047, -1.57075981, -1.03585735,  0.28237642],
        [ 0.37827036,  1.23656925,  2.69465746,  0.00738389, -0.46622797],
        [ 2.61288533,  1.06003315, -1.78241794,  1.14530535, -1.99268424],
        [ 0.89961978, -1.76338331,  0.31624842, -1.53316259,  1.95653327],
        [-0.88055263,  1.93268915,  1.01707394, -0.17647025, -2.35550117]]),
 array([[-0.55940083,  1.97018217, -1.58348281, -1.16246037,  0.16925816],
        [ 0.36988933,  1.06527827,  2.66608288, -0.27695382, -0.72028023],
        [ 2.6142239 ,  1.08739079, -1.77785416,  1.19071819, -1.95210843],
        [ 0.91032536, -1.54458322,  0.35274842, -1.16996123,  2.28104925],
        [-0.86643455,  2.22123367,  1.06520863,  0.30250471, -1.9275429 ]]),
 array([[-0.55940121,  1.97017454, -1.58348408, -1.16247303,  0.16924685],
        [ 0.36988849,  1.06526114,  2.66608002, -0.27698225, -0.72030564],
        [ 2.61422403,  1.08739352, -1.77785371,  1.19072273, -1.95210437],
        [ 0.91032643, -1.54456134,  0.35275207, -1.1699249 ,  2.2810817 ],
        [-0.86643314,  2.22126253,  1.06521344,  0.30255261, -1.9275001 ]]),
 array([[-0.559028  ,  1.97780213, -1.58221165, -1.14981147,  0.18055981],
        [ 0.37072668,  1.08239195,  2.66893777, -0.24854564, -0.69489787],
        [ 2.61409016,  1.08465749, -1.77831013,  1.18618099, -1.95616236],
        [ 0.90925577, -1.56644353,  0.3491017 , -1.20624867,  2.24862686],
        [-0.86784509,  2.19240519,  1.06039949,  0.25465032, -1.97030021]]))

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

x_1 = linalg.inv(A_1.T.dot(A_1)).dot(b)
x_2 = linalg.inv(A_2.T.dot(A_2)).dot(b)
x_3 = linalg.inv(A_3.T.dot(A_3)).dot(b)
x_4 = linalg.inv(A_4.T.dot(A_4)).dot(b)

x_2,3,4 is not very good

x_1, x_2, x_3, x_4

(array([1.3766838 , 2.17441668, 0.71923842, 1.74497268, 2.51336606]),
 array([7.34571152e+06, 1.50131205e+08, 2.50447280e+07, 2.49213142e+08,
        2.22668895e+08]),
 array([-3.50631225e+14, -7.16618086e+15, -1.19545466e+15, -1.18956380e+16,
        -1.06286070e+16]),
 array([  8.64893787, 150.80432815,  25.51352194, 248.46601644,
        222.95560073]))

sigma_plus = linalg.inv(diags([5,4,3,2,1], 0).toarray())
sigma_plus[-1,-1] = 0
sigma_plus

array([[ 0.2       ,  0.        ,  0.        ,  0.        , -0.        ],
       [ 0.        ,  0.25      ,  0.        ,  0.        , -0.        ],
       [ 0.        ,  0.        ,  0.33333333,  0.        , -0.        ],
       [ 0.        ,  0.        ,  0.        ,  0.5       , -0.        ],
       [ 0.        ,  0.        ,  0.        ,  0.        ,  0.        ]])

pseudoinverse

x_pseudo = n.dot(sigma_plus).dot(m).dot(b)

x_pseudo

array([-0.16978402, -0.07097868,  0.99297169, -1.85485295,  1.44494874])
Written on January 4, 2022