Mercurial > matrix-functions

annotate 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
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
0
8f23314345f4 Create local repository for matrix toolboxes. Step #0 done.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff changeset
1 function B = qmult(A)
8f23314345f4 Create local repository for matrix toolboxes. Step #0 done.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff changeset
2 %QMULT Pre-multiply by random orthogonal matrix.
8f23314345f4 Create local repository for matrix toolboxes. Step #0 done.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff changeset
3 % QMULT(A) is Q*A where Q is a random real orthogonal matrix from
8f23314345f4 Create local repository for matrix toolboxes. Step #0 done.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff changeset
4 % the Haar distribution, of dimension the number of rows in A.
8f23314345f4 Create local repository for matrix toolboxes. Step #0 done.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff changeset
5 % Special case: if A is a scalar then QMULT(A) is the same as
8f23314345f4 Create local repository for matrix toolboxes. Step #0 done.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff changeset
6 % QMULT(EYE(A)).
8f23314345f4 Create local repository for matrix toolboxes. Step #0 done.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff changeset
7
8f23314345f4 Create local repository for matrix toolboxes. Step #0 done.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff changeset
8 % Called by RANDSVD.
8f23314345f4 Create local repository for matrix toolboxes. Step #0 done.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff changeset
9 %
8f23314345f4 Create local repository for matrix toolboxes. Step #0 done.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff changeset
10 % Reference:
8f23314345f4 Create local repository for matrix toolboxes. Step #0 done.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff changeset
11 % G.W. Stewart, The efficient generation of random
8f23314345f4 Create local repository for matrix toolboxes. Step #0 done.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff changeset
12 % orthogonal matrices with an application to condition estimators,
8f23314345f4 Create local repository for matrix toolboxes. Step #0 done.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff changeset
13 % SIAM J. Numer. Anal., 17 (1980), 403-409.
8f23314345f4 Create local repository for matrix toolboxes. Step #0 done.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff changeset
14
8f23314345f4 Create local repository for matrix toolboxes. Step #0 done.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff changeset
15 [n, m] = size(A);
8f23314345f4 Create local repository for matrix toolboxes. Step #0 done.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff changeset
16
8f23314345f4 Create local repository for matrix toolboxes. Step #0 done.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff changeset
17 % Handle scalar A.
8f23314345f4 Create local repository for matrix toolboxes. Step #0 done.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff changeset
18 if max(n,m) == 1
8f23314345f4 Create local repository for matrix toolboxes. Step #0 done.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff changeset
19 n = A;
8f23314345f4 Create local repository for matrix toolboxes. Step #0 done.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff changeset
20 A = eye(n);
8f23314345f4 Create local repository for matrix toolboxes. Step #0 done.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff changeset
21 end
8f23314345f4 Create local repository for matrix toolboxes. Step #0 done.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff changeset
22
8f23314345f4 Create local repository for matrix toolboxes. Step #0 done.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff changeset
23 d = zeros(n);
8f23314345f4 Create local repository for matrix toolboxes. Step #0 done.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff changeset
24
8f23314345f4 Create local repository for matrix toolboxes. Step #0 done.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff changeset
25 for k = n-1:-1:1
8f23314345f4 Create local repository for matrix toolboxes. Step #0 done.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff changeset
26
8f23314345f4 Create local repository for matrix toolboxes. Step #0 done.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff changeset
27 % Generate random Householder transformation.
8f23314345f4 Create local repository for matrix toolboxes. Step #0 done.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff changeset
28 x = randn(n-k+1,1);
8f23314345f4 Create local repository for matrix toolboxes. Step #0 done.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff changeset
29 s = norm(x);
8f23314345f4 Create local repository for matrix toolboxes. Step #0 done.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff changeset
30 sgn = sign(x(1)) + (x(1)==0); % Modification for sign(1)=1.
8f23314345f4 Create local repository for matrix toolboxes. Step #0 done.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff changeset
31 s = sgn*s;
8f23314345f4 Create local repository for matrix toolboxes. Step #0 done.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff changeset
32 d(k) = -sgn;
8f23314345f4 Create local repository for matrix toolboxes. Step #0 done.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff changeset
33 x(1) = x(1) + s;
8f23314345f4 Create local repository for matrix toolboxes. Step #0 done.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff changeset
34 beta = s*x(1);
8f23314345f4 Create local repository for matrix toolboxes. Step #0 done.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff changeset
35
8f23314345f4 Create local repository for matrix toolboxes. Step #0 done.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff changeset
36 % Apply the transformation to A.
8f23314345f4 Create local repository for matrix toolboxes. Step #0 done.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff changeset
37 y = x'*A(k:n,:);
8f23314345f4 Create local repository for matrix toolboxes. Step #0 done.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff changeset
38 A(k:n,:) = A(k:n,:) - x*(y/beta);
8f23314345f4 Create local repository for matrix toolboxes. Step #0 done.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff changeset
39
8f23314345f4 Create local repository for matrix toolboxes. Step #0 done.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff changeset
40 end
8f23314345f4 Create local repository for matrix toolboxes. Step #0 done.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff changeset
41
8f23314345f4 Create local repository for matrix toolboxes. Step #0 done.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff changeset
42 % Tidy up signs.
8f23314345f4 Create local repository for matrix toolboxes. Step #0 done.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff changeset
43 for i=1:n-1
8f23314345f4 Create local repository for matrix toolboxes. Step #0 done.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff changeset
44 A(i,:) = d(i)*A(i,:);
8f23314345f4 Create local repository for matrix toolboxes. Step #0 done.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff changeset
45 end
8f23314345f4 Create local repository for matrix toolboxes. Step #0 done.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff changeset
46 A(n,:) = A(n,:)*sign(randn);
8f23314345f4 Create local repository for matrix toolboxes. Step #0 done.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff changeset
47 B = A;