Google Groups no longer supports new Usenet posts or subscriptions. Historical content remains viewable.
Dismiss
rank deficiency tolerance levels in rank() and mldivide()
171 views
Skip to first unread message
Awhan Patnaik
Feb 12, 2012, 5:00:11 AM2/12/12
to
the following piece of code is a simple test case of my actual code
[m, n] = size(A); % in my case m > n
if rank(A) == n
x = A\b; % b is known vector of appropriate dimension
else
disp('rank deficient!!!')
end
i.e. solve the system of equations ONLY when A is full column rank.
i discovered that even when rank() reports a full column rank the backslash operator warns of A being rank deficient!!! on investigating further i gathered that rank() uses svd() to determine the rank while \ uses qr() to solve the over-determined system of equations.
the rank() method uses a tolerance of
tolsvd = max(size(A)) * eps(max(s));
as specified in the corresponding m-file. however there is NO mention of the tolerance used by mldivide or the backslash operator even though it prints a warning when this tolerance level is violated.
i found this link . http://amath.colorado.edu/computing/Matlab/OldTechDocs/ref/qr.html which seems to indicate the tolerance used by \ operator is
tol_mldivide_qr = max(size(A))*eps*abs(R(1,1))
is that right?
could somebody point me to some explanations as to the choice of these bounds? are these standard and published somewhere?
although the two tolerances mentioned above apply to two different objects the rank conclusions about A should be consistent when used on the same computer with same class of variables (double in my case).
also is it possible to pass user specified tolerance to the backslash operator?
in your informed opinion what should be done? should i trust svd() over mldivide() ?
[m, n] = size(A); % in my case m > n
if rank(A) == n
x = A\b; % b is known vector of appropriate dimension
else
disp('rank deficient!!!')
end
i.e. solve the system of equations ONLY when A is full column rank.
i discovered that even when rank() reports a full column rank the backslash operator warns of A being rank deficient!!! on investigating further i gathered that rank() uses svd() to determine the rank while \ uses qr() to solve the over-determined system of equations.
the rank() method uses a tolerance of
tolsvd = max(size(A)) * eps(max(s));
as specified in the corresponding m-file. however there is NO mention of the tolerance used by mldivide or the backslash operator even though it prints a warning when this tolerance level is violated.
i found this link . http://amath.colorado.edu/computing/Matlab/OldTechDocs/ref/qr.html which seems to indicate the tolerance used by \ operator is
tol_mldivide_qr = max(size(A))*eps*abs(R(1,1))
is that right?
could somebody point me to some explanations as to the choice of these bounds? are these standard and published somewhere?
although the two tolerances mentioned above apply to two different objects the rank conclusions about A should be consistent when used on the same computer with same class of variables (double in my case).
also is it possible to pass user specified tolerance to the backslash operator?
in your informed opinion what should be done? should i trust svd() over mldivide() ?
Bruno Luong
Feb 12, 2012, 5:15:10 AM2/12/12
to
"Awhan Patnaik" wrote in message <jh82jb$hht$1...@newscl01ah.mathworks.com>...
Here is a QR base inversion where you can adjust the tolerance:
http://www.mathworks.com/matlabcentral/fileexchange/25453-pseudo-inverse
Bruno
> the following piece of code is a simple test case of my actual code
> in your informed opinion what should be done? should i trust svd() over mldivide() ?
In principle yes, svd method is a bette estimation of rank.
Here is a QR base inversion where you can adjust the tolerance:
http://www.mathworks.com/matlabcentral/fileexchange/25453-pseudo-inverse
Bruno
John D'Errico
Feb 12, 2012, 5:41:10 AM2/12/12
to
"Awhan Patnaik" wrote in message <jh82jb$hht$1...@newscl01ah.mathworks.com>...
> the following piece of code is a simple test case of my actual code
>
> [m, n] = size(A); % in my case m > n
> if rank(A) == n
> x = A\b; % b is known vector of appropriate dimension
> else
> disp('rank deficient!!!')
> end
>
> i.e. solve the system of equations ONLY when A is full column rank.
>
> i discovered that even when rank() reports a full column rank the backslash operator warns of A being rank deficient!!! on investigating further i gathered that rank() uses svd() to determine the rank while \ uses qr() to solve the over-determined system of equations.
>
> the rank() method uses a tolerance of
>
> tolsvd = max(size(A)) * eps(max(s));
>
> as specified in the corresponding m-file. however there is NO mention of the tolerance used by mldivide or the backslash operator even though it prints a warning when this tolerance level is violated.
>
> i found this link . http://amath.colorado.edu/computing/Matlab/OldTechDocs/ref/qr.html which seems to indicate the tolerance used by \ operator is
>
> tol_mldivide_qr = max(size(A))*eps*abs(R(1,1))
>
> is that right?
>
> could somebody point me to some explanations as to the choice of these bounds? are these standard and published somewhere?
>
> although the two tolerances mentioned above apply to two different objects the rank conclusions about A should be consistent when used on the same computer with same class of variables (double in my case).
>
How can a method that does not call svd be ABSOLUTELY
>
> [m, n] = size(A); % in my case m > n
> if rank(A) == n
> x = A\b; % b is known vector of appropriate dimension
> else
> disp('rank deficient!!!')
> end
>
> i.e. solve the system of equations ONLY when A is full column rank.
>
> i discovered that even when rank() reports a full column rank the backslash operator warns of A being rank deficient!!! on investigating further i gathered that rank() uses svd() to determine the rank while \ uses qr() to solve the over-determined system of equations.
>
> the rank() method uses a tolerance of
>
> tolsvd = max(size(A)) * eps(max(s));
>
> as specified in the corresponding m-file. however there is NO mention of the tolerance used by mldivide or the backslash operator even though it prints a warning when this tolerance level is violated.
>
> i found this link . http://amath.colorado.edu/computing/Matlab/OldTechDocs/ref/qr.html which seems to indicate the tolerance used by \ operator is
>
> tol_mldivide_qr = max(size(A))*eps*abs(R(1,1))
>
> is that right?
>
> could somebody point me to some explanations as to the choice of these bounds? are these standard and published somewhere?
>
> although the two tolerances mentioned above apply to two different objects the rank conclusions about A should be consistent when used on the same computer with same class of variables (double in my case).
>
consistent with a method that does call svd?
> also is it possible to pass user specified tolerance to the backslash operator?
>
> in your informed opinion what should be done? should i trust svd() over mldivide() ?
Your results will be highly suspect, and will be terribly
dependent on tiny errors in the least significant bits of
your data.
Yes, trust the signal that rank returns slightly more.
But if backslash is having a problem, that means you
are still at risk of having numerical garbage as a result.
You should consider if something can be improved.
For example, if this is a polynomial model (perhaps the
single largest source of ill-conditioned linear systems)
then you are using too high order polynomials. You
might also consider using centering and scaling of
the problem, but even then, you may still be using
too high an order. If this is from some other scheme,
then you are still treading on thin ice, and should look
for ways to improve the conditioning of your problem.
Just because one method does not flag this system as
ill-posed does not mean it is automatically well-posed!
John
Awhan Patnaik
Feb 14, 2012, 10:04:43 AM2/14/12
to
"Bruno Luong" <b.l...@fogale.findmycountry> wrote in message > > in your informed opinion what should be done? should i trust svd() over mldivide() ?
>
> In principle yes, svd method is a bette estimation of rank.
>
> Here is a QR base inversion where you can adjust the tolerance:
> http://www.mathworks.com/matlabcentral/fileexchange/25453-pseudo-inverse
>
> Bruno
many thanks for your opinion and the link to the code. since i already perform an svd to get a rank estimate i might as well solve my system of equations with svd itself.
> In principle yes, svd method is a bette estimation of rank.
>
> Here is a QR base inversion where you can adjust the tolerance:
> http://www.mathworks.com/matlabcentral/fileexchange/25453-pseudo-inverse
>
> Bruno
Awhan Patnaik
Feb 14, 2012, 10:04:44 AM2/14/12
to
thank you john for the reply.
"John D'Errico" <wood...@rochester.rr.com> wrote in message > > could > >
> >although the two tolerances mentioned above apply to two different objects the rank conclusions about A should be consistent when used on the same computer with same class of variables (double in my case).
> >
> How can a method that does not call svd be ABSOLUTELY
> consistent with a method that does call svd?
>
no john, the expected consistency is with respect to the rank deficiency of the A matrix. i.e. two different methods should be consistent in determining the rank. although i still can not be sure of the tolerance used in mldivide() so assuming the tolerances as in my original post i think that the reason for choosing different tolerances is justified because of the different methods involded. however the rank is a property of the A matrix and inferences about A should be consistent across various methods.
on the other hand if it is known that rank estimates provided by svd methods are currently more trustworthy then either other procedures should use the most trustworthy method or adjust their tolerances to give answers consistent with those of the trustworthy method.
> > in your informed opinion what should be done? should i trust svd() over mldivide() ?
>
> If you are that close to the edge, you should be fearful.
> Your results will be highly suspect, and will be terribly
> dependent on tiny errors in the least significant bits of
> your data.
>
> Yes, trust the signal that rank returns slightly more.
> But if backslash is having a problem, that means you
> are still at risk of having numerical garbage as a result.
>
> You should consider if something can be improved.
> For example, if this is a polynomial model (perhaps the
> single largest source of ill-conditioned linear systems)
> then you are using too high order polynomials. You
> might also consider using centering and scaling of
> the problem, but even then, you may still be using
> too high an order. If this is from some other scheme,
> then you are still treading on thin ice, and should look
> for ways to improve the conditioning of your problem.
>
> Just because one method does not flag this system as
> ill-posed does not mean it is automatically well-posed!
>
> John
i think i should use the condition number of A to guide me on this one. if i examine the ratio of max_singular_value/other_singular_values and truncate after a ratio of about 1000. may be that will work.
another problem that i have is that A is composed in this way
A(i,:) = kron(w(i, :), [1 x(i,:)] %% the 1 in there stems from the use of a bias term
I remember from regression theory that before centering and scaling one removes the column of 1 from the design matrix, however I donot know what I should do in this case.
any ways thanks john for guiding me.
"John D'Errico" <wood...@rochester.rr.com> wrote in message > > could > >
> >although the two tolerances mentioned above apply to two different objects the rank conclusions about A should be consistent when used on the same computer with same class of variables (double in my case).
> >
> How can a method that does not call svd be ABSOLUTELY
> consistent with a method that does call svd?
>
on the other hand if it is known that rank estimates provided by svd methods are currently more trustworthy then either other procedures should use the most trustworthy method or adjust their tolerances to give answers consistent with those of the trustworthy method.
> > in your informed opinion what should be done? should i trust svd() over mldivide() ?
>
> If you are that close to the edge, you should be fearful.
> Your results will be highly suspect, and will be terribly
> dependent on tiny errors in the least significant bits of
> your data.
>
> Yes, trust the signal that rank returns slightly more.
> But if backslash is having a problem, that means you
> are still at risk of having numerical garbage as a result.
>
> You should consider if something can be improved.
> For example, if this is a polynomial model (perhaps the
> single largest source of ill-conditioned linear systems)
> then you are using too high order polynomials. You
> might also consider using centering and scaling of
> the problem, but even then, you may still be using
> too high an order. If this is from some other scheme,
> then you are still treading on thin ice, and should look
> for ways to improve the conditioning of your problem.
>
> Just because one method does not flag this system as
> ill-posed does not mean it is automatically well-posed!
>
> John
another problem that i have is that A is composed in this way
A(i,:) = kron(w(i, :), [1 x(i,:)] %% the 1 in there stems from the use of a bias term
I remember from regression theory that before centering and scaling one removes the column of 1 from the design matrix, however I donot know what I should do in this case.
any ways thanks john for guiding me.
0 new messages