function A = random_positive_definite_matrix(integerN)
  %
  % Generate a random, symmetric positive-definite (SPD) matrix.
  %
  % Since all (SPD) matrices are diagonalizable and have positive
  % eigenvalues, it seems likely that we can generate an SPD matrix by
  % combining random orthogonal matrices with a diagonal matrix of
  % random eigenvalues.
  %
  % I have no proof/evidence that this approach is sound.
  %
  % INPUT:
  %
  %   - ``integerN`` -- The dimension of the resulting matrix.
  %
  % OUTPUT:
  %
  %   - ``A`` -- A symmetric, positive definite matrix.
  %
  U = random_orthogonal_matrix(integerN);
  d = unifrnd(eps, realmax, 1, integerN);
  D = diag(d);
  A = U*D*U';
end
