randsym.m

function S = randsym(n, N)
% Generates random symmetric matrices with normal entries.
% 
% function S = randsym(n)
% function S = randsym(n, N)
%
% S is an n-by-n-by-N array where each slice S(:, :, i) for i = 1..N is a
% random symmetric matrix with upper triangular entries distributed
% independently following a normal distribution (Gaussian, zero mean, unit
% variance).
%
% By default, N = 1.
%
% See also: randrot randskew randherm randskewh

% This file is part of Manopt: www.manopt.org.
% Original author: Nicolas Boumal, Oct. 23, 2018.
% Contributors: 
% Change log: 
%       Oct. 23, 2018 (NB):
%           This is not technically necessary for the rotations factory,
%           but it is counter-intuitive to have access to a function called
%           randskew yet not have one for randsym.
%       June 19, 2019 (NB):
%           Now handles the case n = 1 properly.

    if nargin < 2
        N = 1;
    end
    
    if n == 1
        S = randn(1, 1, N);
        return;
    end

    % Subindices of the (strictly) upper triangular entries of an n-by-n
    % matrix.
    [I, J] = find(triu(ones(n), 1));
    
    K = repmat(1:N, n*(n-1)/2, 1);
    
    % Indices of the strictly upper triangular entries of all N slices of
    % an n-by-n-by-N array.
    L = sub2ind([n n N], repmat(I, N, 1), repmat(J, N, 1), K(:));
    
    % Allocate memory for N random symmetric matrices of size n-by-n and
    % populate each upper triangular entry with a random number following a
    % normal distribution and copy them on the corresponding lower
    % triangular side.
    S = zeros(n, n, N);
    S(L) = randn(size(L));
    S = S + multitransp(S);
    
    % Now populate the diagonal entries:
    
    % Subindices of the diagonal entries of an n-by-n matrix.
    [I, J] = find(eye(n));
    
    K = repmat(1:N, n, 1);
    
    % Indices of the diagonal entries of all N slices of an n-by-n-by-N
    % array.
    L = sub2ind([n n N], repmat(I, N, 1), repmat(J, N, 1), K(:));
    
    S(L) = randn(size(L));
    
end