1function err = test_svd (A) 2%TEST_SVD test factorize(A,'svd') and factorize(A,'cod') for a given matrix 3% 4% Example 5% err = test_svd (A) ; 6% 7% See also test_all 8 9% Copyright 2011-2012, Timothy A. Davis, http://www.suitesparse.com 10 11fprintf ('.') ; 12 13if (nargin < 1) 14 % has rank 3 15 A = magic (4) ; 16end 17 18[m, n] = size (A) ; 19err = 0 ; 20 21for st = 0:1 22 23 if (st == 0) 24 F = factorize (A, 'svd') ; 25 Acond = cond (F) ; 26 else 27 F = factorize (A, 'cod') ; 28 end 29 Apinv = pinv (full (A)) ; 30 31 if (st == 0) 32 assert (ismethod (F, 'norm')) ; 33 assert (ismethod (F, 'pinv')) ; 34 Anorm = norm (full (A)) ; 35 Ainvnorm = norm (pinv (full (A))) ; 36 Fnorm = norm (F) ; 37 e = abs (Anorm - Fnorm) ; 38 Anorm = max (Anorm, Fnorm) ; 39 if (Anorm > 0) 40 e = e / Anorm ; 41 end 42 err = check_err (err, e) ; 43 Anorm = max (Anorm, 1) ; 44 end 45 46 % if B=pinv(A), then A*B*A=A and B*A*B=B 47 B = inverse (F) ; 48 Bnorm = norm (double (B), 1) ; 49 err = check_err (err, ... 50 norm (double(A*B*A) - A, 1) / (Anorm^2 * Bnorm)) ; 51 err = check_err (err, ... 52 norm (double(B*A*B) - double(B), 1) / (Anorm * Bnorm^2)) ; 53 54 for bsparse = 0:1 55 b = rand (m, 1) ; 56 if (bsparse) 57 b = sparse (b) ; 58 end 59 x = Apinv*b ; 60 y = F\b ; 61 if (st == 0) 62 z = pinv(F)*b ; 63 else 64 z = inverse (F)*b ; 65 end 66 if (st == 0 || Acond < 1e13) 67 % skip this for COD for very ill-conditioned problems 68 x = double (x) ; 69 y = double (y) ; 70 z = double (z) ; 71 err = check_err (err, norm (x - y) / (Anorm * norm (x) + norm (b))); 72 err = check_err (err, norm (x - z) / (Anorm * norm (x) + norm (b))); 73 end 74 75 c = rand (1, n) ; 76 if (bsparse) 77 c = sparse (c) ; 78 end 79 x = c*Apinv ; 80 y = c/F ; 81 if (st == 0) 82 z = c*pinv(F) ; 83 else 84 z = c*inverse(F) ; 85 end 86 if (st == 0 || Acond < 1e13) 87 % skip this for COD for very ill-conditioned problems 88 x = double (x) ; 89 y = double (y) ; 90 z = double (z) ; 91 err = check_err (err, norm (x - y) / (Anorm * norm (x) + norm (b))); 92 err = check_err (err, norm (x - z) / (Anorm * norm (x) + norm (b))); 93 end 94 end 95 96 if (st == 0) 97 assert (ismethod (F, 'rank')) ; 98 arank = rank (full (A)) ; 99 if (arank ~= rank (F)) 100 fprintf ('\nrank of A: %d, rank of F: %d\n', arank, rank (F)) ; 101 fprintf ('singular values:\n') ; 102 s1 = svd (A) ; 103 s2 = F.Factors.S ; 104 disp ([s1 s2]) 105 error ('rank mismatch!') ; 106 end 107 108 assert (ismethod (F, 'cond')) ; 109 c1 = cond (full (A)) ; 110 c2 = cond (F) ; 111 if (rank (F) == min (m,n)) 112 e = abs (c1 - c2) ; 113 if (max (c1, c2) > 0) 114 e = e / max (c1, c2) ; 115 end 116 % fprintf ('cond err: %g\n', e) ; 117 err = check_err (err, e) ; 118 else 119 % both c1 and c2 should be large 120 tol = max (m,n) * eps (norm (F)) ; 121 assert (c1 > norm (F) / tol) ; 122 assert (c2 > norm (F) / tol) ; 123 end 124 125 Z = null (F) ; 126 e = norm (A*Z) / Anorm ; 127 % fprintf ('null space err %g (%g)\n', e, err) ; 128 err = check_err (err, e) ; 129 130 [U1, S1, V1] = svd (F) ; 131 [U2, S2, V2] = svd (full (A)) ; 132 err = check_err (err, norm (U1-U2)) ; 133 err = check_err (err, norm (S1-S2) / Anorm) ; 134 err = check_err (err, norm (V1-V2)) ; 135 err = check_err (err, norm (U1*S1*V1'-U2*S2*V2') / Anorm) ; 136 137 [U1, S1, V1] = svd (inverse (F)) ; 138 [U2, S2, V2] = svd (pinv (full (A))) ; 139 err = check_err (err, norm (S1-S2) / Ainvnorm) ; 140 err = check_err (err, norm (U1*S1*V1'-U2*S2*V2') / Ainvnorm) ; 141 142 [U1, S1, V1] = svd (F, 0) ; 143 [U2, S2, V2] = svd (full (A), 0) ; 144 err = check_err (err, norm (U1-U2)) ; 145 err = check_err (err, norm (S1-S2) / Anorm) ; 146 err = check_err (err, norm (V1-V2)) ; 147 err = check_err (err, norm (U1*S1*V1'-U2*S2*V2') / Anorm) ; 148 149 [U1, S1, V1] = svd (F, 'econ') ; 150 [U2, S2, V2] = svd (full (A), 'econ') ; 151 err = check_err (err, norm (U1-U2)) ; 152 err = check_err (err, norm (S1-S2) / Anorm) ; 153 err = check_err (err, norm (V1-V2)) ; 154 err = check_err (err, norm (U1*S1*V1'-U2*S2*V2') / Anorm) ; 155 156 [U1, S1, V1] = svd (F, 'rank') ; 157 err = check_err (err, norm (U1*S1*V1'-U2*S2*V2') / Anorm) ; 158 159 % fprintf ('svd err %g (%g)\n', e, err) ; 160 161 % test the p-norm 162 for k = 0:3 163 if (k == 0) 164 p = 'inf' ; 165 elseif (k == 3) 166 p = 'fro' ; 167 else 168 p = k ; 169 end 170 171 n1 = norm (full (A), p) ; 172 n2 = norm (F, p) ; 173 e = abs (n1 - n2) ; 174 if (n1 > 1) 175 e = e / n1 ; 176 end 177 % fprintf ('norm (A,%d): %g\n', k, e) ; 178 err = check_err (err, e) ; 179 180 n1 = norm (full (A'), p) ; 181 n2 = norm (F', p) ; 182 e = abs (n1 - n2) ; 183 if (n1 > 1) 184 e = e / n1 ; 185 end 186 % fprintf ('norm (A'',%d): %g\n', k, e) ; 187 err = check_err (err, e) ; 188 189 n1 = norm (Apinv, p) ; 190 n2 = norm (pinv (F), p) ; 191 e = abs (n1 - n2) ; 192 if (n1 > 1) 193 e = e / n1 ; 194 end 195 % fprintf ('norm (pinv(A),%d): %g\n', k, e) ; 196 err = check_err (err, e) ; 197 198 n1 = norm (Apinv', p) ; 199 n2 = norm (pinv (F)', p) ; 200 e = abs (n1 - n2) ; 201 if (n1 > 1) 202 e = e / n1 ; 203 end 204 % fprintf ('norm (pinv(A)'',%d): %g\n', k, e) ; 205 err = check_err (err, e) ; 206 end 207 end 208end 209 210function err = check_err (err, e) 211err = max (err, e) ; 212if (err > 1e-6) 213 error ('%g error too high!\n', err) ; 214end 215