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

### Discover what MATLAB® can do for your career.

Opportunities for recent engineering grads.

# 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 ... > 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" wrote in message ... > > > 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 ... > "Florin Neacsu" wrote in message ... > > > > > > 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 ... > 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" wrote in message ... > "Matt J" wrote in message ... > > 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 ... > > > 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 ... > "Bruno Luong" wrote in message ... > > "Matt J" wrote in message ... > > > 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" wrote in message ... > > > 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" wrote in message ... > > > 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 ... > > >> 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" wrote in message ... > "Matt J" wrote in message ... > > > > > >> 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 ... > "Bruno Luong" wrote in message ... > > > > > > 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)
 Subject: Rank deficiency in backslash operations From: Matt J Date: 7 Apr, 2011 01:54:05 Message: 14 of 26 "Florin Neacsu" wrote in message ... > > 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)> 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 ... > "Florin Neacsu" wrote in message ... > > > > 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) ================ > > 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
 Subject: Rank deficiency in backslash operations From: Matt J Date: 7 Apr, 2011 10:40:04 Message: 16 of 26 "Florin Neacsu" wrote in message ... > > > 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
 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 " wrote in message news:inik30\$pou\$1@fred.mathworks.com... > "Florin Neacsu" wrote in message > ... >> >> >> 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" wrote in message ... > 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 ... > "Florin Neacsu" wrote in message ... > > > 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 =================== > > 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 ... > "Matt J" wrote in message ... > > "Florin Neacsu" wrote in message ... > > > > > 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 > =================== > > > > 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" wrote in message news:inl7dn\$3ib\$1@fred.mathworks.com... > "Matt J" wrote in message ... >> "Matt J" wrote in message ... >> > "Florin Neacsu" wrote in message >> > ... >> > >> > > 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 > > =================== >> > >> > 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" wrote in message ... > > > 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" wrote in message ... > 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 " wrote in message news:innv7s\$a3n\$1@fred.mathworks.com... > "Bobby Cheng" wrote in message > ... >> 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" wrote in message ... > 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 " wrote in message news:inpipq\$8t3\$1@fred.mathworks.com... > "Bobby Cheng" wrote in message > ... >> 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.