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