Posted December 6, 2023
The MATLAB builtin function gallery('randsvd',...) can only generate a random matrix or a symmetric positive definite matrix. However, sometimes, we need to test our algorithms for symmetric matrices. Although you can get a symmetric matrix by A = randn(n); A = A + A';, it would be nice to inherit the features from gallery('randsvd'), such as the ability to control the distributions of the singular values. The attached code achieves such desire by
qmult, andWe modified the source code of randsvd, and using two additional functions: qmult and mysign, which are all from Prof. Nicholas J. Higham.
function A = my_randsvd(n, kappa, mode)
%MY_RANDSVD Customized gallery('randsvd')
% A = MY_RANDSVD(N, KAPPA, MODE) generates
% a random symmetric matrix A of order N with
% cond(A) = KAPPA and singular values from distribution MODE.
%
% Input :
% N : Size of the matrix.
%
% KAPPA : pre-defined 2-norm condition number for the output,
% if kappa is negative, it will return a symmetric
% positive definite matrix.
% Default value is 100.
% https://www.mathworks.com/help/matlab/ref/gallery.html
%
%
% MODE : one of the following values:
% 1: one large singular value,
% 2: one small singular value,
% 3: geometrically distributed singular values,
% 4: arithmetically distributed singular values,
% 5: random singular values with uniformly distributed
% logarithm.
% Default value is 3.
%
% Reference:
% The MATLAB 'gallery' function.
classname = 'double';
switch nargin
case 1
kappa = 100;
mode = 3;
case 2
mode = 3;
end
% check if we need export p.d. matrix
if sign(kappa) == 1, pd = false; else, pd = true; end
kappa = abs(kappa);
switch mode % Set up vector of singular values
case 1
sigma = ones(n,1)./kappa;
sigma(1) = 1;
case 2
sigma = ones(n,1);
sigma(n) = 1/kappa;
case 3
factor = kappa^(-1/(n-1));
sigma = factor.^(0:n-1);
case 4
sigma = ones(n,1) - (0:n-1)'/(n-1)*(1-1/kappa);
case 5 % In this case cond(A) <= kappa.
sigma = exp( -rand(n,1)*log(kappa) );
otherwise
error(message('MATLAB:randsvd:invalidMode'));
end
sigma = cast(sigma,classname);
if ~pd % randomly introducing signs
[l,j] = size(sigma(2:n-1));
signs = sign(rand(1,n-2)*2-1)';
signs = reshape(signs,[l,j]);
sigma(2:n-1) = sigma(2:n-1) .* signs;
end
Q = qmult(n,1,classname);
A = Q*diag(sigma)*Q';
A = (A + A')/2; % ensure symmetry
end
Here is the code for qmult and mysign:
function B = qmult(A,method,classname)
%QMULT Pre-multiply matrix by random orthogonal matrix.
% QMULT(A) returns Q*A where Q is a random real orthogonal matrix
% from the Haar distribution of dimension the number of rows in A.
% Special case: if A is a scalar then QMULT(A) is the same as QMULT(EYE(A)).
% QMULT(A,METHOD) specifies how the computations are carried out.
% METHOD = 0 is the default, while METHOD = 1 uses a call to QR,
% which is much faster for large dimensions, even though it uses more flops.
% Called by RANDCOLU, RANDCORR, RANDJORTH, RANDSVD.
% Reference:
% G. W. Stewart, The efficient generation of random
% orthogonal matrices with an application to condition estimators,
% SIAM J. Numer. Anal., 17 (1980), 403-409.
%
% Nicholas J. Higham
% Copyright 1984-2020 The MathWorks, Inc.
if nargin < 2 || isempty(method)
method = 0;
end
% Handle scalar A.
if isscalar(A)
n = A;
A = eye(n,classname);
else
n = size(A, 1);
end
if isempty(A) % nothing to do.
B = A;
return
end
if method == 1
[Q,R] = qr(randn(n));
B = Q*diag(sign(diag(R)))*A;
return
end
d = zeros(n,1,classname);
for k = n-1:-1:1
% Generate random Householder transformation.
x = randn(n-k+1,1);
s = norm(x);
sgn = mysign(x(1));
s = sgn*s;
d(k) = -sgn;
x(1) = x(1) + s;
beta = s*x(1);
% Apply the transformation to A.
y = x'*A(k:n,:);
A(k:n,:) = A(k:n,:) - x*(y/beta);
end
% Tidy up signs.
for i=1:n-1
A(i,:) = d(i)*A(i,:);
end
A(n,:) = A(n,:)*mysign(randn);
B = A;
end
function S = mysign(A)
%MYSIGN True sign function with MYSIGN(0) = 1.
% Called by various matrices in elmat/private.
%
% Nicholas J. Higham
% Copyright 1984-2013 The MathWorks, Inc.
S = sign(A);
S(S==0) = 1;
end