Mercurial > matrix-functions

comparison toolbox/qmult.m @ 2:c124219d7bfa draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Re-add the 1995 toolbox after noticing the statement in the ~higham/mctoolbox/ webpage.
author Antonio Pino Robles <data.script93@gmail.com>
date Thu, 07 May 2015 18:36:24 +0200
parents 8f23314345f4
children
comparison
equal deleted inserted replaced
1:e471a92d17be 2:c124219d7bfa
1 function B = qmult(A)
2 %QMULT Pre-multiply by random orthogonal matrix.
3 % QMULT(A) is Q*A where Q is a random real orthogonal matrix from
4 % the Haar distribution, of dimension the number of rows in A.
5 % Special case: if A is a scalar then QMULT(A) is the same as
6 % QMULT(EYE(A)).
7
8 % Called by RANDSVD.
9 %
10 % Reference:
11 % G.W. Stewart, The efficient generation of random
12 % orthogonal matrices with an application to condition estimators,
13 % SIAM J. Numer. Anal., 17 (1980), 403-409.
14
15 [n, m] = size(A);
16
17 % Handle scalar A.
18 if max(n,m) == 1
19 n = A;
20 A = eye(n);
21 end
22
23 d = zeros(n);
24
25 for k = n-1:-1:1
26
27 % Generate random Householder transformation.
28 x = randn(n-k+1,1);
29 s = norm(x);
30 sgn = sign(x(1)) + (x(1)==0); % Modification for sign(1)=1.
31 s = sgn*s;
32 d(k) = -sgn;
33 x(1) = x(1) + s;
34 beta = s*x(1);
35
36 % Apply the transformation to A.
37 y = x'*A(k:n,:);
38 A(k:n,:) = A(k:n,:) - x*(y/beta);
39
40 end
41
42 % Tidy up signs.
43 for i=1:n-1
44 A(i,:) = d(i)*A(i,:);
45 end
46 A(n,:) = A(n,:)*sign(randn);
47 B = A;