function M = draw_from_gaussian(mu, Sigma, N)
%
% M = draw_from_gaussian(mean, cov, n)
% Each row of M is a random vector drawn from a Gaussian distribution.
% The number of row, N, is 1 if omitted.

if nargin==2
  N = 1;
end

% If Y = CX, Var(Y) = C Var(X) C'.
% So if Var(X)=I, and we want Var(Y)=Sigma, we need to find C. s.t. Sigma = C C'.
% Since Sigma is psd, we have Sigma = U D U' = (U D^0.5) (D'^0.5 U').

mu = mu(:);
n=length(mu);
[U,D,V] = svd(Sigma);
M = randn(n,N);
M = (U*sqrt(D))*M + mu*ones(1,N); % transform each column

M = M';

