matrix pinv( matrix a)		/* pseudo-inverse is  V * pinv(S) * U' */
{
    matrix u, s, v;
    int rank;

    rank = svd( a, u, s, v);

    if( rank == 0)
	return zeros( cols(a), rows(a));
    else
	return v(:,1:rank) * diag( 1/diag(s(1:rank,1:rank))) * u(:,1:rank)';
}
