Matrix Condition and the QR Decomposition - Cleve Moler on Mathematics and Computing

[poster]

Matrix Condition and the QR Decompositionhtml,body,div,span,applet,object,iframe,h1,h2,h3,h4,h5,h6,p,blockquote,pre,a,abbr,acronym,address,big,cite,code,del,dfn,em,font,img,ins,kbd,q,s,samp,small,strike,strong,tt,var,b,u,i,center,dl,dt,dd,ol,ul,li,fieldset,form,label,legend,table,caption,tbody,tfoot,thead,tr,th,td{margin:0;padding:0;border:0;outline:0;font-size:100%;vertical-align:baseline;background:transparent}body{line-height:1}ol,ul{list-style:none}blockquote,q{quotes:none}blockquote:before,blockquote:after,q:before,q:after{content:'';content:none}:focus{outine:0}ins{text-decoration:none}del{text-decoration:line-through}table{border-collapse:collapse;border-spacing:0}html { min-height:100%; margin-bottom:1px; }html body { height:100%; margin:0px; font-family:Arial, Helvetica, sans-serif; font-size:10px; color:#000; line-height:140%; background:#fff none; overflow-y:scroll; }html body td { vertical-align:top; text-align:left; }h1 { padding:0px; margin:0px 0px 25px; font-family:Arial, Helvetica, sans-serif; font-size:1.5em; color:#d55000; line-height:100%; font-weight:normal; }h2 { padding:0px; margin:0px 0px 8px; font-family:Arial, Helvetica, sans-serif; font-size:1.2em; color:#000; font-weight:bold; line-height:140%; border-bottom:1px solid #d6d4d4; display:block; }h3 { padding:0px; margin:0px 0px 5px; font-family:Arial, Helvetica, sans-serif; font-size:1.1em; color:#000; font-weight:bold; line-height:140%; }a { color:#005fce; text-decoration:none; }a:hover { color:#005fce; text-decoration:underline; }a:visited { color:#004aa0; text-decoration:none; }p { padding:0px; margin:0px 0px 20px; }img { padding:0px; margin:0px 0px 20px; border:none; }p img, pre img, tt img, li img, h1 img, h2 img { margin-bottom:0px; }ul { padding:0px; margin:0px 0px 20px 23px; list-style:square; }ul li { padding:0px; margin:0px 0px 7px 0px; }ul li ul { padding:5px 0px 0px; margin:0px 0px 7px 23px; }ul li ol li { list-style:decimal; }ol { padding:0px; margin:0px 0px 20px 0px; list-style:decimal; }ol li { padding:0px; margin:0px 0px 7px 23px; list-style-type:decimal; }ol li ol { padding:5px 0px 0px; margin:0px 0px 7px 0px; }ol li ol li { list-style-type:lower-alpha; }ol li ul { padding-top:7px; }ol li ul li { list-style:square; }.content { font-size:1.2em; line-height:140%; padding: 20px; }pre, code { font-size:12px; }tt { font-size: 1.2em; }pre { margin:0px 0px 20px; }pre.codeinput { padding:10px; border:1px solid #d3d3d3; background:#f7f7f7; }pre.codeoutput { padding:10px 11px; margin:0px 0px 20px; color:#4c4c4c; }pre.error { color:red; }@media print { pre.codeinput, pre.codeoutput { word-wrap:break-word; width:100%; } }span.keyword { color:#0000FF }span.comment { color:#228B22 }span.string { color:#A020F0 }span.untermstring { color:#B20000 }span.syscmd { color:#B28C00 }span.typesection { color:#A0522D }.footer { width:auto; padding:10px 0px; margin:25px 0px 0px; border-top:1px dotted #878787; font-size:0.8em; line-height:140%; font-style:italic; color:#878787; text-align:left; float:none; }.footer p { margin:0px; }.footer a { color:#878787; }.footer a:hover { color:#878787; text-decoration:underline; }.footer a:visited { color:#878787; }table th { padding:7px 5px; text-align:left; vertical-align:middle; border: 1px solid #d6d4d4; font-weight:bold; }table td { padding:7px 5px; text-align:left; vertical-align:top; border:1px solid #d6d4d4; } The QR decomposition provides an estimate of the matrix condition number.

Contents

• Query
• Short answer
• QR Decomposition
• tricond
• qr function
• Gallery
• Tests
• Plot

Query

Professor Magdy Hanna of the Department of Engineering Mathematics and Physics at Fayoum University in Fayoum, Egypt, asked me:

If we have the QR decomposition of a matrix, with or without pivoting, is there an easy way to get the condition number of the matrix without starting from scratch and using the SVD?Short answer

The short answer to Prof. Hanna's question is:

Yes. The ratio of the largest and smallest diagonal elements of qr(A) is often a good estimate of cond(A).Compare this estimate to the exact condition of a magic square. Produce a score measuring how well the estimate does.

A = magic(11); R = qr(A); d = abs(diag(R)); est = max(d)/min(d); kappa = cond(A); disp(" est kappa score") disp([est kappa est/kappa]) est kappa score 2.6964 11.1021 0.2429QR Decomposition

A more complete answer involves column pivoting.

The QR decomposition without column pivoting is an orthogonal Q and an upper triangular R so that

A = Q*RThe QR composition with column pivoting is an orthogonal Q, an upper triangular R and a permutation P that reorders the columns of A so that

A*P = Q*RIn either case, Q and P are orthogonal, so

cond(R) = cond(A)Any estimate of the condition of R provides an estimate of the condition of A.

tricond

An estimate the condition of a triangular matrix is the ratio of the largest and smallest elements on the diagonal.

type tricond function cond = tricond(T) % Estimate condition of lower or upper triangular matrix. if istril(T) || istriu(T) D = diag(abs(T)); cond = max(D)/min(D); else cond = NaN; end endqr function

Here is the 4-by-4 Pascal matrix

A = pascal(4)A = 1 1 1 1 1 2 3 4 1 3 6 10 1 4 10 20This matrix is reasonably well conditioned.

kappa = cond(A)kappa = 691.9374With one output, the qr function provides the "Q-less" QR decomposition with no pivoting.

(With two outputs, the qr function produces the same R and also reveals the orthogonal Q.)

R = qr(A)R = -2.0000 -5.0000 -10.0000 -17.5000 0 -2.2361 -6.7082 -14.0872 0 0 1.0000 3.5000 0 0 0 -0.2236For this matrix tricond without pivoting produces a poor estimate of the exact condition.

trico = tricond(R); kappa = cond(A); disp(" tricond kappa score") disp([trico kappa trico/kappa]) tricond kappa score 10.0000 691.9374 0.0145With three outputs, qr does column pivoting and returns a different R.

[Q,R,P] = qr(A); RR = -22.7376 -5.2336 -1.5393 -12.0065 0 -1.6153 -1.2034 -1.3387 0 0 0.4270 -0.2170 0 0 0 -0.0638Compare tricond and the exact condition.

trico = tricond(R); disp(" tricond kappa score") disp([trico kappa trico/kappa]) tricond kappa score 356.6259 691.9374 0.5154Column pivoting produces a much better tricond estimate for this matrix.

Gallery

Let's do more tests using these matrices from the MATLAB Gallery.

funs = ["moler", "parter", "binomial", "chebspec", "cauchy",... "chebspec", "chebvand","circul", "frank", "ris"];One of these Gallery matrices is the "moler" matrix, even though I didn't invent it. A precursor to the "moler" matrix is a badly conditioned upper triangular matrix with ones on the diagonal, minus ones above the diagonal and zeros below.

n = 5; U = eye(n) - triu(ones(n),1)U = 1 -1 -1 -1 -1 0 1 -1 -1 -1 0 0 1 -1 -1 0 0 0 1 -1 0 0 0 0 1The "moler" matrix is

M = U'*UM = 1 -1 -1 -1 -1 -1 2 0 0 0 -1 0 3 1 1 -1 0 1 4 2 -1 0 1 2 5The "moler" matrix squares the condition of U. The condition of M grows like 4^n.

kappa = cond(M)kappa = 865.9801The tricond score without column pivoting is very poor.

R = qr(M); trico = tricond(R); disp(" tricond kappa score") disp([trico kappa trico/kappa]) tricond kappa score 20.7364 865.9801 0.0239With pivoting the score is much better.

[~,R,~] = qr(M); trico = tricond(R); disp(" tricond kappa score") disp([trico kappa trico/kappa]) tricond kappa score 555.0198 865.9801 0.6409These scores are in the first line of the following output and in the upper left-hand corners of the 3-D bar graphs.

Tests

[Dnopiv,Dpiv,xt,yt] = tester(funs); fun n nopiv piv moler 5 0.024 0.641 10 0.000 0.426 15 0.000 0.337 20 0.000 0.287 25 0.000 0.372 parter 5 0.610 0.804 10 0.564 0.776 15 0.536 0.757 20 0.517 0.744 25 0.502 0.735 binomial 5 0.795 0.926 10 0.387 0.817 15 0.150 0.768 20 0.056 0.696 25 0.021 0.674 chebspec 5 0.120 0.110 10 0.236 0.095 15 0.183 0.358 20 0.132 0.165 25 0.022 0.318 cauchy 5 0.217 0.376 10 0.034 0.242 15 0.016 0.247 20 0.058 0.116 25 0.102 0.222 chebspec 5 0.120 0.110 10 0.236 0.095 15 0.183 0.358 20 0.132 0.165 25 0.022 0.318 chebvand 5 0.055 0.577 10 0.002 0.326 15 0.000 0.338 20 0.000 0.492 25 0.000 0.082 circul 5 0.370 0.370 10 0.252 0.252 15 0.209 0.209 20 0.180 0.180 25 0.162 0.162 frank 5 0.390 0.390 10 0.257 0.296 15 0.191 0.248 20 Inf 0.013 25 0.039 0.100 ris 5 0.610 0.804 10 0.564 0.776 15 0.536 0.757 20 0.517 0.744 25 0.502 0.735 Plot

[c,ax1,ax2] = ploter(funs,Dnopiv,Dpiv,xt,yt);
Published with MATLAB® R2026a





More...