A MATLAB function to generate random symmetric matrix
Published:
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
- first, generate singular values,
- then construct a random orthogonal matrix using
qmult, and - finally multiply the orthogonal matrix at both sides of the diagonal matrix that contains these singular values.
We modified the source code of randsvd, and using two additional functions: qmult and mysign, which are all from Prof. Nicholas J. Higham.
- Version: Thursday December 6, 2023
- Acknowledgement: 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