Mercurial > matrix-functions

comparison toolbox/randsvd.m @ 0:8f23314345f4 draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Create local repository for matrix toolboxes. Step #0 done.
author Antonio Pino Robles <data.script93@gmail.com>
date Wed, 06 May 2015 14:56:53 +0200
parents
children
comparison
equal deleted inserted replaced
-1:000000000000 0:8f23314345f4
1 function A = randsvd(n, kappa, mode, kl, ku)
2 %RANDSVD Random matrix with pre-assigned singular values.
3 % RANDSVD(N, KAPPA, MODE, KL, KU) is a (banded) random matrix of order N
4 % with COND(A) = KAPPA and singular values from the distribution MODE.
5 % N may be a 2-vector, in which case the matrix is N(1)-by-N(2).
6 % Available types:
7 % MODE = 1: one large singular value,
8 % MODE = 2: one small singular value,
9 % MODE = 3: geometrically distributed singular values,
10 % MODE = 4: arithmetically distributed singular values,
11 % MODE = 5: random singular values with unif. dist. logarithm.
12 % If omitted, MODE defaults to 3, and KAPPA defaults to SQRT(1/EPS).
13 % If MODE < 0 then the effect is as for ABS(MODE) except that in the
14 % original matrix of singular values the order of the diagonal entries
15 % is reversed: small to large instead of large to small.
16 % KL and KU are the lower and upper bandwidths respectively; if they
17 % are omitted a full matrix is produced.
18 % If only KL is present, KU defaults to KL.
19 % Special case: if KAPPA < 0 then a random full symmetric positive
20 % definite matrix is produced with COND(A) = -KAPPA and
21 % eigenvalues distributed according to MODE.
22 % KL and KU, if present, are ignored.
23
24 % Reference:
25 % N.J. Higham, Accuracy and Stability of Numerical Algorithms,
26 % Society for Industrial and Applied Mathematics, Philadelphia, PA,
27 % USA, 1996; sec. 26.3.
28
29 % This routine is similar to the more comprehensive Fortran routine xLATMS
30 % in the following reference:
31 % J.W. Demmel and A. McKenney, A test matrix generation suite,
32 % LAPACK Working Note #9, Courant Institute of Mathematical Sciences,
33 % New York, 1989.
34
35 if nargin < 2, kappa = sqrt(1/eps); end
36 if nargin < 3, mode = 3; end
37 if nargin < 4, kl = n-1; end % Full matrix.
38 if nargin < 5, ku = kl; end % Same upper and lower bandwidths.
39
40 if abs(kappa) < 1, error('Condition number must be at least 1!'), end
41 posdef = 0; if kappa < 0, posdef = 1; kappa = -kappa; end % Special case.
42
43 p = min(n);
44 m = n(1); % Parameter n specifies dimension: m-by-n.
45 n = n(max(size(n)));
46
47 if p == 1 % Handle case where A is a vector.
48 A = randn(m, n);
49 A = A/norm(A);
50 return
51 end
52
53 j = abs(mode);
54
55 % Set up vector sigma of singular values.
56 if j == 3
57 factor = kappa^(-1/(p-1));
58 sigma = factor.^[0:p-1];
59
60 elseif j == 4
61 sigma = ones(p,1) - (0:p-1)'/(p-1)*(1-1/kappa);
62
63 elseif j == 5 % In this case cond(A) <= kappa.
64 sigma = exp( -rand(p,1)*log(kappa) );
65
66 elseif j == 2
67 sigma = ones(p,1);
68 sigma(p) = 1/kappa;
69
70 elseif j == 1
71 sigma = ones(p,1)./kappa;
72 sigma(1) = 1;
73 end
74
75 % Convert to diagonal matrix of singular values.
76 if mode < 0
77 sigma = sigma(p:-1:1);
78 end
79 sigma = diag(sigma);
80
81 if posdef % Handle special case.
82 Q = qmult(p);
83 A = Q'*sigma*Q;
84 A = (A + A')/2; % Ensure matrix is symmetric.
85 return
86 end
87
88 if m ~= n
89 sigma(m, n) = 0; % Expand to m-by-n diagonal matrix.
90 end
91
92 if kl == 0 & ku == 0 % Diagonal matrix requested - nothing more to do.
93 A = sigma;
94 return
95 end
96
97 % A = U*sigma*V, where U, V are random orthogonal matrices from the
98 % Haar distribution.
99 A = qmult(sigma');
100 A = qmult(A');
101
102 if kl < n-1 | ku < n-1 % Bandwidth reduction.
103 A = bandred(A, kl, ku);
104 end