www.gusucode.com > elmat工具箱matlab源码程序 > elmat/private/randcolu.m

    function A = randcolu(x, m, k, classname)
%RANDCOLU Random matrix with normalized columns and specified singular values.
%   A = GALLERY('RANDCOLU',N) is a random N-by-N matrix with columns of
%   unit 2-norm, with random singular values whose squares are from a
%   uniform distribution.  GALLERY('RANDCOLU',N,M), where M >= N,
%   produces an M-by-N matrix.
%   A'*A is a correlation matrix of the form produced by
%   GALLERY('RANDCORR',N).
%
%   GALLERY('RANDCOLU',X), where X is an N-vector (N > 1), produces
%   a random N-by-N matrix having singular values given by the vector X.
%   X must have nonnegative elements whose sum of squares is N.
%   GALLERY('RANDCOLU',X,M), where M >= N, produces an M-by-N matrix.
%   GALLERY('RANDCOLU',X,M,K) provides a further option:
%   For K = 0 (the default) DIAG(X) is initially subjected to a random
%       two-sided orthogonal transformation and then a sequence of
%       Givens rotations is applied.
%   For K = 1, the initial transformation is omitted. This is much faster,
%       but the resulting matrix may have zero entries.

%  Reference:
%  P. I. Davies and N. J. Higham, Numerically stable generation of
%  correlation matrices and their factors, BIT, 40 (2000), pp. 640-651.
%
%   Nicholas J. Higham
%   Copyright 1984-2005 The MathWorks, Inc.

if isempty(k), k = 0; end
if length(x) == 1
   n = x;
   x = rand(n,1);
   x = sqrt(n)*x/norm(x);
   if isempty(m), m = n; end
   x = cast(x,classname);
else
   n = length(x);
   if isempty(m), m = n; end
end

if m < n
  error(message('MATLAB:randcolu:SizeM'))
end
if abs(sum(x.^2) - n)/n > 100*eps(classname) | any(x < 0)
   error(message('MATLAB:randcolu:InvalidX'))
end

A = zeros(m,n,classname);
for i = 1:n
   A(i,i) = x(i);
end

if k == 0
   % Forming A --> U*A*V where U and V are two random orthogonal matrices.
   A = qmult(A,[],classname);
   A = qmult(A',[],classname)';
end

a = sum(A .* A);

y = find(a<1);
z = find(a>1);

while length(y) > 0 && length(z) > 0

   i = y(ceil(rand*length(y)));
   j = z(ceil(rand*length(z)));
   if i > j, temp = i; i = j; j = temp; end

   aij = A(:,i)'*A(:,j);
   alpha = sqrt(aij^2 - (a(i)-1)*(a(j)-1));

   t(1) = (aij + mysign(aij)*alpha)/(a(j)-1);
   t(2) = (a(i)-1)/((a(j)-1)*t(1));
   t = t(ceil(rand*2));  % Choose randomly from the two roots.
   c = 1/sqrt(1 + t^2) ;  s = t*c ;

   A(:,[i,j]) =   A(:,[i,j]) * [c s ; -s c] ;

   a(j) = a(j) + a(i) - 1;
   a(i) = 1;

   y = find(a<1);
   z = find(a>1);

end