How does the MATLAB A\B compute, when A is singular square and full (dense) matrix?

user_6650451 0
Viewed 95

I am working on linear equations like A*X = B. I got a matrix like A = magic(4) and B = [1;3;2;4]. I used several methods to solve this problem, also including MATLAB A\B.

What's annoyed is the resulf of MATLAB A\B is different from those of NumPy, SciPy and LAPACK (LAPACKE_dgetrf, LAPACKE_dgels, LAPACKE_dgelsd).

The folowing are the results.

MATLAB A\B (Ubuntu Matlab 2021a and Matlab Online)

>> A = magic(4)


A =

    16     2     3    13
     5    11    10     8
     9     7     6    12
     4    14    15     1

>> B = [1;3;2;4]

B =

     1
     3
     2
     4

>> A \ B
RCOND =  4.625929e-18。  

ans =

   -0.0098
    0.0735
    0.1985
    0.0319
>> lsqr(A, B)

ans =

    0.0110
    0.1360
    0.1360
    0.0110

>> x = pinv(A)*B

x =

    0.0110
    0.1360
    0.1360
    0.0110

>> [L, U, P] = lu(A) %lu decomposition

L =

    1.0000         0         0         0
    0.2500    1.0000         0         0
    0.5625    0.4352    1.0000         0
    0.3125    0.7685    1.0000    1.0000


U =

   16.0000    2.0000    3.0000   13.0000
         0   13.5000   14.2500   -2.2500
         0         0   -1.8889    5.6667
         0         0         0    0.0000


P =

     1     0     0     0
     0     0     0     1
     0     0     1     0
     0     1     0     0
>> y = L \ (P * B) % LU tril solver

y =

    1.0000
    3.7500
   -0.1944
    0.0000

>> x = U \ y % LU triu solver

x =

   -0.0098
    0.0735
    0.1985
    0.0319

Python numpy.linalg.solve(a, b)

import numpy as np

a = np.array([[16., 2., 3., 13.], [5., 11., 10., 8.], [9., 7., 6., 12.], [4., 14., 15., 1.]]);
b = np.array([[1.],[3.],[2.],[4]]);
x = np.linalg.solve(a, b)
'''
x = 
array([[ 0.10799632],
       [ 0.42693015],
       [-0.15487132],
       [-0.0859375 ]])
'''

Python scipy.linalg.solve(a, b)

from scipy import linalg

x1 = linalg.solve(a, b)
'''
x1 = 
array([[ 0.10799632],
       [ 0.42693015],
       [-0.15487132],
       [-0.0859375 ]])
'''

Octave Online

octave:1> A = magic(4)
A =

   16    2    3   13
    5   11   10    8
    9    7    6   12
    4   14   15    1

octave:2> B = [1;3;2;4]
B =

   1
   3
   2
   4

octave:3> A\B
warning: matrix singular to machine precision, rcond = 1.30614e-17
ans =

   0.011029
   0.136029
   0.136029
   0.011029

Lapack LAPACKE_degelsd

a = 
  Matrix(Row = 4, Col = 4, Major = ColMajor)
            16             2             3            13
             5            11            10             8
             9             7             6            12
             4            14            15             1

B = 
  Vector(Size = 4)
             1
             3
             2
             4

info = 0
X = 
  Vector(Size = 4)
      0.107996
       0.42693
     -0.154871
    -0.0859375

My problem is what the MATLAB A \ B exactly computes. I know that pinv(A) * B computes the minimum norm least-squares solution, and so do np.linalg.solve and scipy.linalg.solve.

And the MATLAB reference page said:

  • If A is a square matrix, then A\B is roughly equal to inv(A)*B, but MATLAB processes A\B differently and more robustly.

  • If the rank of A is less than the number of columns in A, then x = A\B is not necessarily the minimum norm solution. You can compute the minimum norm least-squares solution using x = lsqminnorm(A,B) or x = pinv(A)*B.

Well I know the solution is not unique when A is nearly singular. There are many feasible solutions (like numpy, octave and LAPACK. They solve this problem to minimize the norm of ||AX - B|| by QR, SVD(pinv)). However what makes matlab get its solution, and what's the magic?

So, I want to ask what the operation A\B exactly computes and what the result of A\B is when A is a singular square matrix.

P.S.

I do know there isn't an unique solution for linear equations when its coefficient matrix is singular or nearly singular .aka. ill conditioned. According to the reference page of \(mldivide), x = A\B is not necessarily the minimum norm solution. I want to know the underlying details of A\B. If it is not necessarily the minimum norm solution, then what? As I know, scipy, numpy, LAPACK and Octave get a minimum norm solution. Is there another fomulation for this problem that I haven't know?

linear-algebra matrix matlab numpy python
1 Answers

There are several ways to try to see this problem, but the one I would recommend is that it makes no sense to try to solve Ax = B, when A is ill-conditioned:

Look at the 1D case for example. What is x when 0.x = b?

We know that this problem has infinitely many solutions if b = 0 or none if b~=0. What does MATLAB say about that:

0\0

ans =

   NaN


0\1

ans =

   Inf

Not very convincing. What about lsqminnorm:

lsqminnorm(0,0)

ans =

     0

lsqminnorm(0,1)

ans =

     0

See? lsqminnorm gives you a valid answer only when there are infinitely many solutions. When there are none, it just tries to minimize || 0.x - b ||. Does this make sense? Are you satisfied with it returning x = 0 as a solution of 0.x = 1?

Back to your problem. Let's try a different value for B:

B2 = [1;1;2;1];

out2 = A\B2;

A*out2

out_minnorm = lsqminnorm(A,B2);

A*out_minnorm

out_pinv = pinv(A)*B2;

A*out_pinv

ans =

         0
   -1.0000
         0
   -2.6250


ans =

    1.1500
    1.4500
    1.5500
    0.8500


ans =

    1.1500
    1.4500
    1.5500
    0.8500

Are any of these values for x acceptable? Even though the last two results are the minimum norm solutions, these don't really mean much.

So, my take on your problem is to check if rcond gets close to machine precision, because not returning any result is in my opinion better than returning something that doesn't make sense:

if rcond(A) < 1e-12

    error('A is singular, the problem seems to be ill-conditioned');

else

    out = B\A;

end
user_43830320
Related