Got Questions? Get Answers.
Discover MakerZone

MATLAB and Simulink resources for Arduino, LEGO, and Raspberry Pi

Learn more

Discover what MATLAB® can do for your career.

Opportunities for recent engineering grads.

Apply Today

Thread Subject:
Rank deficiency in backslash operations

Subject: Rank deficiency in backslash operations

From: Matt J

Date: 5 Apr, 2011 23:00:19

Message: 1 of 26

Alright, time to broaden, hopefully, my understanding of numerical linear algebra.

I have a 1899960 x 6 matrix AA, with cond(AA)= 21.5308
which seems fairly well-conditioned. I also have a data vector, Y, of type single.

When I do a backslash operation, I get a rank deficiency warning

K>> x=AA\Y
Warning: Rank deficient, rank = 5, tol = 2.2649e-001.

x =

  1.0e+003 *

    2.1030
         0
   -0.0030
   -0.0095
   -0.0095
   -0.0095


When I convert Y to double type and repeat, the warning goes away and I get significantly different values

K>> x=AA\double(Y)

x =

  1.0e+003 *

    2.9346
   -0.8393
    0.0006
   -0.0060
   -0.0060
   -0.0060


So my questions are

(1) Don't the rank warnings refer to the rank of A? If so, why does Y being single or double make a difference? Surely the precision of Y shouldn't affect the rank of A....

(2) Why do I have rank issues at all when AA is so well-conditioned?

Subject: Rank deficiency in backslash operations

From: Florin Neacsu

Date: 5 Apr, 2011 23:13:04

Message: 2 of 26

"Matt J" wrote in message <ing6u2$o0p$1@fred.mathworks.com>...
> Alright, time to broaden, hopefully, my understanding of numerical linear algebra.
>
> I have a 1899960 x 6 matrix AA, with cond(AA)= 21.5308
> which seems fairly well-conditioned. I also have a data vector, Y, of type single.
>
> When I do a backslash operation, I get a rank deficiency warning
>
> K>> x=AA\Y
> Warning: Rank deficient, rank = 5, tol = 2.2649e-001.
>
> x =
>
> 1.0e+003 *
>
> 2.1030
> 0
> -0.0030
> -0.0095
> -0.0095
> -0.0095
>
>
> When I convert Y to double type and repeat, the warning goes away and I get significantly different values
>
> K>> x=AA\double(Y)
>
> x =
>
> 1.0e+003 *
>
> 2.9346
> -0.8393
> 0.0006
> -0.0060
> -0.0060
> -0.0060
>
>
> So my questions are
>
> (1) Don't the rank warnings refer to the rank of A? If so, why does Y being single or double make a difference? Surely the precision of Y shouldn't affect the rank of A....
>
> (2) Why do I have rank issues at all when AA is so well-conditioned?

Hello,

From the doc of mldivide :

" Data Type Support
When computing X = A\B or X = A/B,
the matrices A and B can have data type double or single.
The following rules determine the data type of the result:If both A and B have
type double, X has type double.If either A or B has
type single, X has type single. "

In your case when Y is single, AA becomes single too, hence the issues with the rank of AA.

Regards,
Florin

Subject: Rank deficiency in backslash operations

From: Matt J

Date: 6 Apr, 2011 05:29:05

Message: 3 of 26

"Florin Neacsu" <fneacsu2@gmail.com> wrote in message <ing7m0$5th$1@fred.mathworks.com>...
>
>
> From the doc of mldivide :
>
> " Data Type Support
> When computing X = A\B or X = A/B,
> the matrices A and B can have data type double or single.
> The following rules determine the data type of the result:If both A and B have
> type double, X has type double.If either A or B has
> type single, X has type single. "
>
> In your case when Y is single, AA becomes single too, hence the issues with the rank of AA.
=================

No, that doesn't quite explain it to. me. The doc says that the result X will be cast to type single. That doesn't seem to explain if/why the input operands like AA are converted to singles as well...

Even if it's true that AA becomes single, it still doesn't why AA is rank deficient when the condition number is so small. When I convert AA to single explicitly, the condition number is barely affected:

 K>> cond(AA)

ans =

   21.5308

K>> cond(single(AA))

ans =

   21.5304

Subject: Rank deficiency in backslash operations

From: Florin Neacsu

Date: 6 Apr, 2011 16:58:04

Message: 4 of 26

"Matt J" wrote in message <ingtn1$p9t$1@fred.mathworks.com>...
> "Florin Neacsu" <fneacsu2@gmail.com> wrote in message <ing7m0$5th$1@fred.mathworks.com>...
> >
> >
> > From the doc of mldivide :
> >
> > " Data Type Support
> > When computing X = A\B or X = A/B,
> > the matrices A and B can have data type double or single.
> > The following rules determine the data type of the result:If both A and B have
> > type double, X has type double.If either A or B has
> > type single, X has type single. "
> >
> > In your case when Y is single, AA becomes single too, hence the issues with the rank of AA.
> =================
>
> No, that doesn't quite explain it to. me. The doc says that the result X will be cast to type single. That doesn't seem to explain if/why the input operands like AA are converted to singles as well...
>
> Even if it's true that AA becomes single, it still doesn't why AA is rank deficient when the condition number is so small. When I convert AA to single explicitly, the condition number is barely affected:
>
> K>> cond(AA)
>
> ans =
>
> 21.5308
>
> K>> cond(single(AA))
>
> ans =



Hi,

If you look whar rank does, you will notice that the tolerence is computed with respect to the size of your matrix and the data type.

When you do AA\Y you get

Warning: Rank deficient, rank = 5, tol = 2.2649e-001.

If you compute
>1899960*eps('single')

ans =

    0.2265
This implies that AA is considered to be single here. I guess Y being single, matlab treats AA as single also.
But if you try

>1899960*eps('double')

ans =

  4.2188e-010
so the tolerence is a lot smaller which will increase the rank.

Which would explain why in the second case AA\double(Y) you don't get the warning.

Regards,
Florin

Subject: Rank deficiency in backslash operations

From: Bruno Luong

Date: 6 Apr, 2011 18:08:02

Message: 5 of 26

"Matt J" wrote in message <ing6u2$o0p$1@fred.mathworks.com>...
> Alright, time to broaden, hopefully, my understanding of numerical linear algebra.
>
> I have a 1899960 x 6 matrix AA, with cond(AA)= 21.5308
> which seems fairly well-conditioned. I also have a data vector, Y, of type single.

Odd, what do you get when you run
rank(AA)

5?

It does seem like a bug to me.

Bruno

Subject: Rank deficiency in backslash operations

From: Matt J

Date: 6 Apr, 2011 18:45:38

Message: 6 of 26

"Bruno Luong" <b.luong@fogale.findmycountry> wrote in message <inia62$3bo$1@fred.mathworks.com>...
> "Matt J" wrote in message <ing6u2$o0p$1@fred.mathworks.com>...
> > Alright, time to broaden, hopefully, my understanding of numerical linear algebra.
> >
> > I have a 1899960 x 6 matrix AA, with cond(AA)= 21.5308
> > which seems fairly well-conditioned. I also have a data vector, Y, of type single.
>
> Odd, what do you get when you run
> rank(AA)
>
> 5?
=========================

Not with AA, but with single(AA), yes, which is consistent with Florin's remarks.

>> rank(AA)

ans =

     6

>> rank(single(AA))

ans =

     5

This alone wouldn't bother me so much. It might mean only that the tolerance used in RANK computations is quite stringent.

What bothers me more is the inconsistency with cond(AA)=21.5308
I know the following to be the true solution

xtrue =

  1.0e+003 *

    2.9346
   -0.8393
    0.0006
   -0.0060
   -0.0060
   -0.0060

There is a relative error in the data of


>> DataError = norm(AA*xtrue-Y)/norm(AA*xtrue)

DataError =

  4.6729e-008

There is a relative error in the solution of

>> xsol=AA\Y; SolutionError=norm(xsol-xtrue)/norm(xtrue)
Warning: Rank deficient, rank = 5, tol = 2.2649e-001.

SolutionError =

    0.3871


As I understand the theory, cond(AA) should provide a bound on

>> SolutionError/DataError

ans =

  8.2839e+006


but of course this is much greater than 21.53.

Subject: Rank deficiency in backslash operations

From: Bruno Luong

Date: 6 Apr, 2011 20:05:20

Message: 7 of 26

"Matt J" wrote in message <inicch$bd9$1@fred.mathworks.com>...
>
>
> What bothers me more is the inconsistency with cond(AA)=21.5308
> I know the following to be the true solution
>

Agree, a (single) matrix that does not have a full rank must not have condition of 21.

There is a bug in cond.

The rest seems to be consistent.

Bruno

Subject: Rank deficiency in backslash operations

From: Florin Neacsu

Date: 6 Apr, 2011 20:28:04

Message: 8 of 26

"Matt J" wrote in message <inicch$bd9$1@fred.mathworks.com>...
> "Bruno Luong" <b.luong@fogale.findmycountry> wrote in message <inia62$3bo$1@fred.mathworks.com>...
> > "Matt J" wrote in message <ing6u2$o0p$1@fred.mathworks.com>...
> > > Alright, time to broaden, hopefully, my understanding of numerical linear algebra.
> > >
> > > I have a 1899960 x 6 matrix AA, with cond(AA)= 21.5308
> > > which seems fairly well-conditioned. I also have a data vector, Y, of type single.
> >
> > Odd, what do you get when you run
> > rank(AA)
> >
> > 5?
> =========================
>
> Not with AA, but with single(AA), yes, which is consistent with Florin's remarks.
>
> >> rank(AA)
>
> ans =
>
> 6
>
> >> rank(single(AA))
>
> ans =
>
> 5
>
> This alone wouldn't bother me so much. It might mean only that the tolerance used in RANK computations is quite stringent.
>
> What bothers me more is the inconsistency with cond(AA)=21.5308
> I know the following to be the true solution
>
> xtrue =
>
> 1.0e+003 *
>
> 2.9346
> -0.8393
> 0.0006
> -0.0060
> -0.0060
> -0.0060
>
> There is a relative error in the data of
>
>
> >> DataError = norm(AA*xtrue-Y)/norm(AA*xtrue)
>
> DataError =
>
> 4.6729e-008
>
> There is a relative error in the solution of
>
> >> xsol=AA\Y; SolutionError=norm(xsol-xtrue)/norm(xtrue)
> Warning: Rank deficient, rank = 5, tol = 2.2649e-001.
>
> SolutionError =
>
> 0.3871
>
>
> As I understand the theory, cond(AA) should provide a bound on
>
> >> SolutionError/DataError
>
> ans =
>
> 8.2839e+006
>
>
> but of course this is much greater than 21.53.

Hi,

Could you tell us what cond(AA,inf) outputs ? cond(AA) uses the 2 norm and I've only seen proofs for the bound you are using in the case of inf norm. It makes sense that such a bound be true for all norms applyiable to matrices, but ...

Florin

Subject: Rank deficiency in backslash operations

From: Matt J

Date: 6 Apr, 2011 20:57:04

Message: 9 of 26

"Florin Neacsu" <fneacsu2@gmail.com> wrote in message <iniick$qlp$1@fred.mathworks.com>...
>
>
> Could you tell us what cond(AA,inf) outputs ?
===============

Nothing of value. cond(A,inf) only supports square A

>> cond(AA,inf)
??? Error using ==> cond at 36
A is rectangular. Use the 2 norm.




>cond(AA) uses the 2 norm and I've only seen proofs for the bound you are using in the case of inf norm. It makes sense that such a bound be true for all norms applyiable to matrices, but ...
=====================

It's true for all induced norms. See here for a generic derivation

http://en.wikipedia.org/wiki/Condition_number#Matrices

Subject: Rank deficiency in backslash operations

From: Matt J

Date: 6 Apr, 2011 21:08:05

Message: 10 of 26

"Bruno Luong" <b.luong@fogale.findmycountry> wrote in message <inih20$4c2$1@fred.mathworks.com>...
>
>
> Agree, a (single) matrix that does not have a full rank must not have condition of 21.
>
> There is a bug in cond.
===========

Then there is a bug in SVD.

>> s=svd(single(AA))

s =

    1.5383
    0.9909
    0.9909
    0.9745
    0.8456
    0.0714

>> max(s)/min(s)

ans =

   21.5304




I think the tol used by RANK is just way too high. How is the following threshold appropriate for singular values ranging from .0714 to 1.5383?

>> tol = max(size(single(AA))) * eps(norm(single(AA)))

tol =

    0.2265

Subject: Rank deficiency in backslash operations

From: Bruno Luong

Date: 6 Apr, 2011 21:24:05

Message: 11 of 26

"Matt J" wrote in message <iniknl$7f5$1@fred.mathworks.com>...

>
> >> tol = max(size(single(AA))) * eps(norm(single(AA)))

For a while, I also wonder the term max(size(single(AA))) in estimate the rank.

Here is my guess: The relative large error you get is mainly due to the sum of a long vectors.

Bruno

Subject: Rank deficiency in backslash operations

From: Matt J

Date: 6 Apr, 2011 22:25:04

Message: 12 of 26

"Bruno Luong" <b.luong@fogale.findmycountry> wrote in message <inilll$o6a$1@fred.mathworks.com>...
> "Matt J" wrote in message <iniknl$7f5$1@fred.mathworks.com>...
>
> >
> > >> tol = max(size(single(AA))) * eps(norm(single(AA)))
>
> For a while, I also wonder the term max(size(single(AA))) in estimate the rank.
>
> Here is my guess: The relative large error you get is mainly due to the sum of a long vectors.
=====================

That certainly seems to have an influence, as does AA being converted internally to single, but the most major factor seems to be the rank test done by MLDIVIDE. I've done several experiments below using a QR approach (which bypasses the rank test) so that I could study the effect of each independently.

Even though not all of the ErrorRatios below obey the conditional number bound, they all do much better than MLDIVIDE.
 


>> [Q,R]=qr(single(AA),0); Q=double(Q); R=double(R);
>> xsol=R\(Q'*double(Y)); ErrorRatio=norm(double(xsol)-xtrue)/norm(xtrue)/DataError

ErrorRatio =

   44.3266




>> [Q,R]=qr(AA,0);
>> xsol=R\(Q'*Y); ErrorRatio=norm(double(xsol)-xtrue)/norm(xtrue)/DataError

ErrorRatio =

   77.7574



>> [Q,R]=qr(AA,0);
>> xsol=R\(Q'*double(Y)); ErrorRatio=norm(double(xsol)-xtrue)/norm(xtrue)/DataError

ErrorRatio =

    0.0060

Subject: Rank deficiency in backslash operations

From: Florin Neacsu

Date: 6 Apr, 2011 22:40:22

Message: 13 of 26

"Matt J" wrote in message <iniknl$7f5$1@fred.mathworks.com>...
> "Bruno Luong" <b.luong@fogale.findmycountry> wrote in message <inih20$4c2$1@fred.mathworks.com>...
> >
> >
> > Agree, a (single) matrix that does not have a full rank must not have condition of 21.
> >
> > There is a bug in cond.
> ===========
>
> Then there is a bug in SVD.
>
> >> s=svd(single(AA))
>
> s =
>
> 1.5383
> 0.9909
> 0.9909
> 0.9745
> 0.8456
> 0.0714
>
> >> max(s)/min(s)
>
> ans =
>
> 21.5304
>
>
>
>
> I think the tol used by RANK is just way too high. How is the following threshold appropriate for singular values ranging from .0714 to 1.5383?
>
> >> tol = max(size(single(AA))) * eps(norm(single(AA)))
>
> tol =
>
> 0.2265

Hi,

The last value of s < tol which yields rank = 5. Is then AA considered singular or not ? If it is singular, can cond(A) still equal max(s)/min(s) ? Is min(s) = 0.0714 or 0, since min(s)<tol ?

@Matt J : Sorry for asking you what cond(AA,inf) yields. It was stupid of me. It is confusing how this thing works and I got lost in the details and forgot that AA is not square.

Florin

Subject: Rank deficiency in backslash operations

From: Matt J

Date: 7 Apr, 2011 01:54:05

Message: 14 of 26

"Florin Neacsu" <fneacsu2@gmail.com> wrote in message <iniq4m$6go$1@fred.mathworks.com>...
>
> The last value of s < tol which yields rank = 5. Is then AA considered singular or not ? If it is singular, can cond(A) still equal max(s)/min(s) ? Is min(s) = 0.0714 or 0, since min(s)<tol ?
================

I wouldn't rely on this rank estimate to assess whether AA is to be considered singular or not. Let's take a simpler example,

>> Q=diag([1, 2^-4]), Qs=single(Q);

Q =

    1.0000 0
         0 0.0625


These matrices consist of zeros and small powers of two. They therefore represent exactly the quantities intended here with no floating point error, either in double precision. One consequence of this is that MATLAB considers them equal:


>> isequal(Q,Qs)

ans =

     1

Now let's pad these matrices with lots of rows of zeros

>> Q(1899960,2)=0; Qs(1899960,2)=0;

Surely you'd agree that such a padding operation should not change the rank of these matrices when no floating point error has been introduced. And...since MATLAB still considers these matrices equal,

>> isequal(Q,Qs)

ans =

     1

surely it should also continue to at least think that their ranks are equal. It doesn't however:


>> rank(Q), rank(Qs)

ans =

     2


ans =

     1

Subject: Rank deficiency in backslash operations

From: Florin Neacsu

Date: 7 Apr, 2011 02:20:23

Message: 15 of 26

"Matt J" wrote in message <inj5ft$1ag$1@fred.mathworks.com>...
> "Florin Neacsu" <fneacsu2@gmail.com> wrote in message <iniq4m$6go$1@fred.mathworks.com>...
> >
> > The last value of s < tol which yields rank = 5. Is then AA considered singular or not ? If it is singular, can cond(A) still equal max(s)/min(s) ? Is min(s) = 0.0714 or 0, since min(s)<tol ?
> ================
>
> I wouldn't rely on this rank estimate to assess whether AA is to be considered singular or not. Let's take a simpler example,
>
> >> Q=diag([1, 2^-4]), Qs=single(Q);
>
> Q =
>
> 1.0000 0
> 0 0.0625
>
>
> These matrices consist of zeros and small powers of two. They therefore represent exactly the quantities intended here with no floating point error, either in double precision. One consequence of this is that MATLAB considers them equal:
>
>
> >> isequal(Q,Qs)
>
> ans =
>
> 1
>
> Now let's pad these matrices with lots of rows of zeros
>
> >> Q(1899960,2)=0; Qs(1899960,2)=0;
>
> Surely you'd agree that such a padding operation should not change the rank of these matrices when no floating point error has been introduced. And...since MATLAB still considers these matrices equal,
>
> >> isequal(Q,Qs)
>
> ans =
>
> 1
>
> surely it should also continue to at least think that their ranks are equal. It doesn't however:
>
>
> >> rank(Q), rank(Qs)
>
> ans =
>
> 2
>
>
> ans =
>
> 1

I agree with what you are saying. Once again, this happens because of the tolerance difference in single and double.

My remark was directed towards cond. Shouldn't cond take into account this also, and give : cond(single(AA))=inf ? (AA of original post). That is because one singular value is zero, as it is <tol ...

Regards,
Florin

Subject: Rank deficiency in backslash operations

From: Matt J

Date: 7 Apr, 2011 10:40:04

Message: 16 of 26

"Florin Neacsu" <fneacsu2@gmail.com> wrote in message <inj717$oni$1@fred.mathworks.com>...
>
>
> I agree with what you are saying. Once again, this happens because of the tolerance difference in single and double.
>
> My remark was directed towards cond. Shouldn't cond take into account this also, and give : cond(single(AA))=inf ? (AA of original post). That is because one singular value is zero, as it is <tol ...
===================

I don't think so. What I'm trying to say is that the tol is wrong, and therefore it, and not COND, is suspect (in my mind at least).

Or to put it another way, I don't see a reason why the rank tol would be set up to consider Q rank 1 here simply because I've added rows of zeros. In a world of ideal, perfect precision math, adding rows of zeros should never have this effect. And since this particular Q is perfectly representable in binary floating point, I don't see why ideas from the perfect precision math world wouldn't apply here.

Subject: Rank deficiency in backslash operations

From: Bobby Cheng

Date: 7 Apr, 2011 15:53:44

Message: 17 of 26

Matt,

Since MLDIVIDE is trying to give you a least square solution, how much does
that affect the residual norm norm(A*x-b)?

This is the most relevant measure in accessing the goodness of the computed
solution.

---Bob.


"Matt J " <mattjacREMOVE@THISieee.spam> wrote in message
news:inik30$pou$1@fred.mathworks.com...
> "Florin Neacsu" <fneacsu2@gmail.com> wrote in message
> <iniick$qlp$1@fred.mathworks.com>...
>>
>>
>> Could you tell us what cond(AA,inf) outputs ?
> ===============
>
> Nothing of value. cond(A,inf) only supports square A
>
>>> cond(AA,inf)
> ??? Error using ==> cond at 36
> A is rectangular. Use the 2 norm.
>
>
>
>
>>cond(AA) uses the 2 norm and I've only seen proofs for the bound you are
>>using in the case of inf norm. It makes sense that such a bound be true
>>for all norms applyiable to matrices, but ...
> =====================
>
> It's true for all induced norms. See here for a generic derivation
>
> http://en.wikipedia.org/wiki/Condition_number#Matrices

Subject: Rank deficiency in backslash operations

From: Matt J

Date: 7 Apr, 2011 17:51:05

Message: 18 of 26

"Bobby Cheng" <bcheng@mathworks.com> wrote in message <inkmm8$4m1$1@fred.mathworks.com>...
> Matt,
>
> Since MLDIVIDE is trying to give you a least square solution, how much does
> that affect the residual norm norm(A*x-b)?
>
> This is the most relevant measure in accessing the goodness of the computed
> solution.
==============

Bob, the residual norm obtained with MLDIVIDE is unquestionably poorer. You don't need my data to reproduce this comparison. Below, I provide a test script giving a much simpler example,


%Simulated data
AA=[diag([1, 2^-4]) ; zeros(1899960 ,2)];
Xtrue=[1;pi];
Y=single(AA*Xtrue);
[Q,R]=qr(AA,0);


Xmld=AA\Y; %Solution using MLDIVIDE

Xqr=R\(Q'*Y); %Solution using QR, no rank checking


%%%%%Output


>> Xmld

Xmld =

     1
     0

>> norm(AA*Xmld-Y)

ans =

    0.1963


>> Xqr

Xqr =

    1.0000
    3.1416

>> norm(AA*Xqr-Y)

ans =

     0

Subject: Rank deficiency in backslash operations

From: Matt J

Date: 7 Apr, 2011 19:52:05

Message: 19 of 26

"Matt J" wrote in message <ink4a4$d4f$1@fred.mathworks.com>...
> "Florin Neacsu" <fneacsu2@gmail.com> wrote in message <inj717$oni$1@fred.mathworks.com>...
>
> > My remark was directed towards cond. Shouldn't cond take into account this also, and give : cond(single(AA))=inf ? (AA of original post). That is because one singular value is zero, as it is <tol ...
> ===================
>
> I don't think so. What I'm trying to say is that the tol is wrong, and therefore it, and not COND, is suspect (in my mind at least).
================

Maybe more relevant to your remark is: I might support COND being designed to be consistent with the output of RANK, but that's if I could trust RANK. Below are modified versions of RANK and COND that I might trust:



function [r,s] = rank(A,tol)
%RANK Customized matrix rank function, by Matt J.


s = svd(A);
normA=max(s);
if nargin==1
   tol = norm(A) * eps(normA);
end
r = sum(s > tol);





function c = cond(A, p)
%COND as designed by TMW slightly tweaked by Matt J

if nargin == 1, p = 2; end

if isempty(A)
    c = zeros(class(A));
    return
end
if issparse(A)
    warning('MATLAB:cond:SparseNotSupported', ...
        'Using CONDEST instead of COND for sparse matrix.')
    c = condest(A);
    return
end

[m, n] = size(A);
if m ~= n && p ~= 2
   error('MATLAB:cond:normMismatchSizeA', 'A is rectangular. Use the 2 norm.')
end

if p == 2

  %%%%%%%%%%Matt J changed the next 2 lines

   [r,s] = svd(A);
   if ~r % Handle singular matrix, as determined by RANK ()


       c = Inf(class(A));
   else
       c = max(s)./min(s);
   end
else
% We'll let NORM pick up any invalid p argument.
   c = norm(A, p) * norm(inv(A), p);
end

Subject: Rank deficiency in backslash operations

From: Florin Neacsu

Date: 7 Apr, 2011 20:39:19

Message: 20 of 26

"Matt J" wrote in message <inl4l5$bt2$1@fred.mathworks.com>...
> "Matt J" wrote in message <ink4a4$d4f$1@fred.mathworks.com>...
> > "Florin Neacsu" <fneacsu2@gmail.com> wrote in message <inj717$oni$1@fred.mathworks.com>...
> >
> > > My remark was directed towards cond. Shouldn't cond take into account this also, and give : cond(single(AA))=inf ? (AA of original post). That is because one singular value is zero, as it is <tol ...
> > ===================
> >
> > I don't think so. What I'm trying to say is that the tol is wrong, and therefore it, and not COND, is suspect (in my mind at least).
> ================
>
> Maybe more relevant to your remark is: I might support COND being designed to be consistent with the output of RANK, but that's if I could trust RANK. Below are modified versions of RANK and COND that I might trust:
>
>
>
> function [r,s] = rank(A,tol)
> %RANK Customized matrix rank function, by Matt J.
>
>
> s = svd(A);
> normA=max(s);
> if nargin==1
> tol = norm(A) * eps(normA);
> end
> r = sum(s > tol);
>
>
>
>
>
> function c = cond(A, p)
> %COND as designed by TMW slightly tweaked by Matt J
>
> if nargin == 1, p = 2; end
>
> if isempty(A)
> c = zeros(class(A));
> return
> end
> if issparse(A)
> warning('MATLAB:cond:SparseNotSupported', ...
> 'Using CONDEST instead of COND for sparse matrix.')
> c = condest(A);
> return
> end
>
> [m, n] = size(A);
> if m ~= n && p ~= 2
> error('MATLAB:cond:normMismatchSizeA', 'A is rectangular. Use the 2 norm.')
> end
>
> if p == 2
>
> %%%%%%%%%%Matt J changed the next 2 lines
>
> [r,s] = svd(A);
> if ~r % Handle singular matrix, as determined by RANK ()
>
>
> c = Inf(class(A));
> else
> c = max(s)./min(s);
> end
> else
> % We'll let NORM pick up any invalid p argument.
> c = norm(A, p) * norm(inv(A), p);
> end

Hi,

Regarding rank : The rank of m by n matrix is <=min(m,n). So I wonder why matlab is using
>tol = max(size(A))*eps(max(s)); in their rank function. As you pointed out adding rows of zeros to a matrix, changes the rank when it shouldn't. So you think it might be a typo and that line should actually be
>tol = min(size(A))*eps(max(s)); ? This way, adding lines of zeros won't change a thing.

If I may, in your function normA and norm(A) isn't the same thing ? maybe you could change the line
>tol = norm(A) * eps(normA);
into
>tol = normA * eps(normA);
You avoid calling svd which might be time consuming.
Also, why did you choose norm(A) ? (I don't have any suggestions, I am just curious how you thought through this).
Thanks.
Florin

Subject: Rank deficiency in backslash operations

From: Bobby Cheng

Date: 7 Apr, 2011 22:39:43

Message: 21 of 26

Matt,

I think you have a point and there is an opportunity in improving MATLAB.
For numerical rank, this is always going to be an issue for any particular
tolerence to miss what the users ideally want. For now, you can cast the
single LHS to double first. I don't think of it as a bug but it is something
that we can do better.

I will make a record of this as something to work on and think about. If you
like, you can also contact our support and lay out the case personally.

Regards,
---Bob.

"Florin Neacsu" <fneacsu2@gmail.com> wrote in message
news:inl7dn$3ib$1@fred.mathworks.com...
> "Matt J" wrote in message <inl4l5$bt2$1@fred.mathworks.com>...
>> "Matt J" wrote in message <ink4a4$d4f$1@fred.mathworks.com>...
>> > "Florin Neacsu" <fneacsu2@gmail.com> wrote in message
>> > <inj717$oni$1@fred.mathworks.com>...
>> >
>> > > My remark was directed towards cond. Shouldn't cond take into account
>> > > this also, and give : cond(single(AA))=inf ? (AA of original post).
>> > > That is because one singular value is zero, as it is <tol ...
>> > ===================
>> >
>> > I don't think so. What I'm trying to say is that the tol is wrong, and
>> > therefore it, and not COND, is suspect (in my mind at least).
>> ================
>>
>> Maybe more relevant to your remark is: I might support COND being
>> designed to be consistent with the output of RANK, but that's if I could
>> trust RANK. Below are modified versions of RANK and COND that I might
>> trust:
>>
>>
>>
>> function [r,s] = rank(A,tol)
>> %RANK Customized matrix rank function, by Matt J.
>>
>>
>> s = svd(A);
>> normA=max(s);
>> if nargin==1
>> tol = norm(A) * eps(normA);
>> end
>> r = sum(s > tol);
>>
>>
>>
>>
>>
>> function c = cond(A, p)
>> %COND as designed by TMW slightly tweaked by Matt J
>>
>> if nargin == 1, p = 2; end
>>
>> if isempty(A)
>> c = zeros(class(A));
>> return
>> end
>> if issparse(A)
>> warning('MATLAB:cond:SparseNotSupported', ...
>> 'Using CONDEST instead of COND for sparse matrix.')
>> c = condest(A);
>> return
>> end
>>
>> [m, n] = size(A);
>> if m ~= n && p ~= 2
>> error('MATLAB:cond:normMismatchSizeA', 'A is rectangular. Use the 2
>> norm.')
>> end
>>
>> if p == 2
>>
>> %%%%%%%%%%Matt J changed the next 2 lines
>>
>> [r,s] = svd(A);
>> if ~r % Handle singular matrix, as determined by RANK ()
>>
>>
>> c = Inf(class(A));
>> else
>> c = max(s)./min(s);
>> end
>> else
>> % We'll let NORM pick up any invalid p argument.
>> c = norm(A, p) * norm(inv(A), p);
>> end
>
> Hi,
>
> Regarding rank : The rank of m by n matrix is <=min(m,n). So I wonder why
> matlab is using
>>tol = max(size(A))*eps(max(s)); in their rank function. As you pointed out
>>adding rows of zeros to a matrix, changes the rank when it shouldn't. So
>>you think it might be a typo and that line should actually be
>>tol = min(size(A))*eps(max(s)); ? This way, adding lines of zeros won't
>>change a thing.
>
> If I may, in your function normA and norm(A) isn't the same thing ? maybe
> you could change the line
>>tol = norm(A) * eps(normA);
> into
>>tol = normA * eps(normA);
> You avoid calling svd which might be time consuming. Also, why did you
> choose norm(A) ? (I don't have any suggestions, I am just curious how you
> thought through this).
> Thanks.
> Florin

Subject: Rank deficiency in backslash operations

From: Matt J

Date: 8 Apr, 2011 20:19:05

Message: 22 of 26

"Florin Neacsu" <fneacsu2@gmail.com> wrote in message <inl7dn$3ib$1@fred.mathworks.com>...
>
>
> Regarding rank : The rank of m by n matrix is <=min(m,n). So I wonder why matlab is using
> >tol = max(size(A))*eps(max(s)); in their rank function. =========================

I think I'm beginning to understand. One can argue that an approximate null vector x should satisfy

norm(A*x,p)<=tolerance

The computation on the left hand side, however, can in general suffer a pile up of round-off errors in some semi-proportion to numel(A)*eps(A). The tolerance value on the right hand side would need to account for that somehow. So, it's quite conceivable that TMW reached a size-dependent tolerance value with reasoning along these lines...

These considerations, however, don't take the potential sparsity of A into account. For sparse A, like in my examples, the numerical-noise pile-up can be much lower, and size(A)-dependent error bounds will be much too pessimistic.

 
> If I may, in your function normA and norm(A) isn't the same thing ? maybe you could change the line
> >tol = norm(A) * eps(normA);
> into
> >tol = normA * eps(normA);
> You avoid calling svd which might be time consuming.
===============

Quite right. Thanks!



> Also, why did you choose norm(A) ? (I don't have any suggestions, I am just curious how you thought through this).
=================

I'm no longer as confident as I was in this choice. I basically considered only diagonal A. For a diagonal matrix, the nullity is the number of zero-valued diagonals. For an approximate nullity, it stands to reason you'd take the diagonal values that are small compared to the others. Hence it made sense to threshold relative to max(diag(A)).

Subject: Rank deficiency in backslash operations

From: Matt J

Date: 8 Apr, 2011 21:38:04

Message: 23 of 26

"Bobby Cheng" <bcheng@mathworks.com> wrote in message <inleff$2ck$1@fred.mathworks.com>...
> Matt,
>
> I think you have a point and there is an opportunity in improving MATLAB.
> For numerical rank, this is always going to be an issue for any particular
> tolerence to miss what the users ideally want. For now, you can cast the
> single LHS to double first. I don't think of it as a bug but it is something
> that we can do better.
======================

Thanks, Bob. Yes, I agree that an all-encompassing tolerance could be hard to find.
If my guesses in Message #22 of the reasoning behind the tolerance formula is correct, however, maybe it would be worth considering a modification of it to account for sparsity. Something along the lines, maybe of

tol= max(size(A)) * eps(norm(A)) *nnz(A)/numel(A)

Subject: Rank deficiency in backslash operations

From: Bobby Cheng

Date: 8 Apr, 2011 23:43:12

Message: 24 of 26

Matt,

Does you matrix sparse enough to store in the sparse format?

---Bob.

"Matt J " <mattjacREMOVE@THISieee.spam> wrote in message
news:innv7s$a3n$1@fred.mathworks.com...
> "Bobby Cheng" <bcheng@mathworks.com> wrote in message
> <inleff$2ck$1@fred.mathworks.com>...
>> Matt,
>>
>> I think you have a point and there is an opportunity in improving MATLAB.
>> For numerical rank, this is always going to be an issue for any
>> particular tolerence to miss what the users ideally want. For now, you
>> can cast the single LHS to double first. I don't think of it as a bug but
>> it is something that we can do better.
> ======================
>
> Thanks, Bob. Yes, I agree that an all-encompassing tolerance could be hard
> to find.
> If my guesses in Message #22 of the reasoning behind the tolerance formula
> is correct, however, maybe it would be worth considering a modification of
> it to account for sparsity. Something along the lines, maybe of
>
> tol= max(size(A)) * eps(norm(A)) *nnz(A)/numel(A)

Subject: Rank deficiency in backslash operations

From: Matt J

Date: 9 Apr, 2011 12:18:02

Message: 25 of 26

"Bobby Cheng" <bcheng@mathworks.com> wrote in message <ino6if$2q3$1@fred.mathworks.com>...
> Matt,
>
> Does you matrix sparse enough to store in the sparse format?
>
> ---Bob.
======================

Converting to sparse wouldn't help if the goal is to study rank since RANK doesn't support sparse type. It also doesn't improve AA\single(Y) because MLDIVIDE doesn't support single precision Y operands.

Anyway, I'm not sure that it matters from TMW's perspective. I think you want your RANK function to be able to know that the columns of the identity matrix are linearly independent:

>> Q=eye(1899960*5,2,'single');

>> rank(Q)

ans =

     0

Subject: Rank deficiency in backslash operations

From: Bobby Cheng

Date: 9 Apr, 2011 12:38:47

Message: 26 of 26

I would agree with you. Food for thought, here.

---Bob.

"Matt J " <mattjacREMOVE@THISieee.spam> wrote in message
news:inpipq$8t3$1@fred.mathworks.com...
> "Bobby Cheng" <bcheng@mathworks.com> wrote in message
> <ino6if$2q3$1@fred.mathworks.com>...
>> Matt,
>>
>> Does you matrix sparse enough to store in the sparse format?
>>
>> ---Bob.
> ======================
>
> Converting to sparse wouldn't help if the goal is to study rank since RANK
> doesn't support sparse type. It also doesn't improve AA\single(Y) because
> MLDIVIDE doesn't support single precision Y operands.
>
> Anyway, I'm not sure that it matters from TMW's perspective. I think you
> want your RANK function to be able to know that the columns of the
> identity matrix are linearly independent:
>
>>> Q=eye(1899960*5,2,'single');
>
>>> rank(Q)
>
> ans =
>
> 0
>

Tags for this Thread

What are tags?

A tag is like a keyword or category label associated with each thread. Tags make it easier for you to find threads of interest.

Anyone can tag a thread. Tags are public and visible to everyone.

Contact us