% sampleRandomSings(n,k) calls testRandomSings(n) k times and find the
% frequency with which LU decomposition generates a zero vector in the
% last row and the frequency with which the rank is measured at n-1.
function [p,q] = sampleRandomSings(n,k)
sz = 0;
sr = 0;
for i=1:k
[z,r] = testRandomSing(n);
sz = sz+z;
sr = sr+r;
end
p = sz/k;
q = sr/k;
end
% Construct a random singular n x n matrix by making the nth column a random
% weighted sum of the others. Check whether row echelon decomposition
% gives an exact zero vector in the final row. Test whether rank(u) gives
% the right rank (n-1)
function [t,r] = testRandomSing(n)
a = rand(n,n-1);
b = rand(n-1,1);
a(:,n) = a*b;
[l,u] = lu(a); % u is the row-echelon decomposition of a;
t=(sum(u(n,:)==0)==n); % check that the last row of u is the zero vector
r=(rank(u)==n-1);
end