null() fails for a symbolic numeric 5x5 matrix

Define the following symbolic 5x5 matrix whose entries are all numbers:

a = str2sym('-(sqrt(sqrt(5) + 5)*(69*sqrt(2) - 13*sqrt(10)))/1320');
b = str2sym('-(sqrt(sqrt(5) + 5)*(29*sqrt(2) + 21*sqrt(10)))/1320');
c = str2sym('(sqrt(2)*(4*sqrt(5) + 49)*sqrt(sqrt(5) + 5))/330');
A = [c a b b a; a c a b b; b a c a b; b b a c a; a b b a c];

Then find the null space:

null(A) returns Empty sym: 5-by-0, but null(eval(A)) returns [1 1 1 1 1]' which is the correct result.

As a check, simplify(A*[1 1 1 1 1]') returns [0 0 0 0 0]'.

Ale

Answer 1

[Alessandro Russo]

Just to inform you that in MATLAB version 9.13.0.2049777 (R2022b) the bug seems to have been solved.

Ciao.

Ale

>> a = str2sym('-(sqrt(sqrt(5) + 5)*(69*sqrt(2) - 13*sqrt(10)))/1320');
>> b = str2sym('-(sqrt(sqrt(5) + 5)*(29*sqrt(2) + 21*sqrt(10)))/1320');
>> c = str2sym('(sqrt(2)*(4*sqrt(5) + 49)*sqrt(sqrt(5) + 5))/330');
>> A = [c a b b a; a c a b b; b a c a b; b b a c a; a b b a c];
>> simplify(null(A))
 
ans =
 
1
1
1
1
1
 
>> 

Answer 2

[Bruno Luong]

Just a wild guess but when you call null(A) with symbolic matrix, MATLAB tries to solve for lambda the equation

det(A - lambda*eye(5)) == 0

This is a fifth order polynomial on lambda and MATLAB does not know how to solve it as Galois has teaches us.

[Alessandro Russo]

If you try to compute the eigenvalues of the matrix, the answer is

root(cumbersome 5th degree polynomial in z, z, 1)
...
root(cumbersome 5th degree polynomial in z, z, 5)

so MATLAB realizes that he cannot find the roots of that polynomial. But for null() there are no warnings/hints whatsoever.

Alessandro Russo

[Bruno Luong]

Yeah, I don't believe that my answer makes much sense to be honest.

[Alessandro Russo]

My guess is the following: det(A) gives

>> det(A)
 
ans =
 
((5^(1/2) + 5)^(5/2)*(419541*2^(1/2)*5^(1/2) - 39928*5^(1/2)*10^(1/2) + 199640*2^(1/2) - 419541*10^(1/2)))/9783848250

but after simplification

>> simplify(ans)
 
ans =
 
0

so maybe null() does not think hard enough!!

Ale

[John D'Errico]

No. null is not using a determinant to test to see if the matrix is singular. I would bet a reasonable sum of money on that assertion. (A determinant is almost NEVER a good choice to do much of anything, unless one is assigning it as a homework problem for a student, where determinants seem to be a clear favorite.) Instead, my best bet is null uses rref, a common way one might compute a nullspace.

Answer 3

[John D'Errico]

It is often difficult to know what a compiled code has done under the hood, especially when symbolic operations are involved. However, my guess is that null may be employing rref to compute the nullspace. There are many ways to compute a nullspace, not only using svd as MATLAB does for double arrays.

rref(A)
ans =
[1, 0, 0, 0, 0]
[0, 1, 0, 0, 0]
[0, 0, 1, 0, 0]
[0, 0, 0, 1, 0]
[0, 0, 0, 0, 1]

That would suggest that A has full rank, and so the nullspace is empty. rref appears to be incorrect here.

How about QR?

[Q,R] = qr(A);

R(5,:)
ans =
[0, 0, 0, 0, 0]

So QR is able to see that A has less than full rank, and we could compute a null space from a QR decomposition. For example working on the transpose of A to compute the nullspace of the columns of A using a QR, we would see that

[Q,R] = qr(A.');

simplify(Q(:,end))
ans =
5^(1/2)/5
5^(1/2)/5
5^(1/2)/5
5^(1/2)/5
5^(1/2)/5

But SVD fails, because it runs into the standard 5th degree polynomial problems.

simplify(svd(A))
ans =
root(z^5 + z^4*((2*2^(1/2)*10^(...

Of course, then rank fails, though it too dives under the mex cloak of darkness to do its work.

rank(A)
ans =
5

Is A truly rank 4? It looks to be rank 4, at least in double form.

svd(double(A))
ans =
0.951056516295154
0.951056516295153
0.718891283158874
0.718891283158874
4.84853835847603e-17

Even vpa agrees.

vpa(svd(A))
ans =
7.1125683591344274549613591959346e-38i
0.71889128315887390879918908979534 - 0.000000000000000024094219551455689173676197254085i
0.71889128315887496851465541352151 + 0.000000000000000024094219551455843032602081297669i
0.95105651629515317160385237328668 - 0.000000000000000018212507998940379872651231332919i
0.95105651629515397262902629347034 + 0.0000000000000000182125079989402213862389448737i

Anyway, my guess is that sym/null, which then dives into mupadmex code for the nullspace, probably uses an rref implementation. Hard to know for sure though.

No comment is given from the code. It thinks A is full rank. No warning would be deemed necessary.

I'd just suggest sending this in as a bug report. Note that Answers is NOT official technical support.

Answer 4

[Alessandro Russo]

Hi everybody,

I filed a bug report (Case Number 05453152). Incoherent results also appear with qr() and simplify:

>> [Q,R]=qr(A);
>> simplify(A-Q*R)
 
ans =
 
[0, 0, 0, 0, 0]
[0, 0, 0, 0, 0]
[0, 0, 0, 0, 0]
[0, 0, 0, 0, 0]
[0, 0, 0, 0, 0]

% so the decomposition is correct
 
>> null(A)
 
ans =
 
Empty sym: 5-by-0

% this was the original wrong behaviour

% but..

>> simplify(null(Q*R))
 
ans =
 
1
1
1
1
1

% this time the answer is correct!!

% however:

>> simplify(null(simplify(Q*R)))
 
ans =
 
Empty sym: 5-by-0
 
>>

% it seems that "simplify" misleads "null"!

Have a nice Sunday

Alessandro