JuliaLang/julia — pull request #8859 by zgimbutas (2014-10-30 22:30:34+00:00)
Title: Fix least squares

The least squares routines lose accuracy in Julia 0.3. This branch
applies a set of simple fixes to address this issue. We test the implementation
of least squares by forming an ill-conditioned operator A, a right hand
side b in the range of A, and finding a least squares solution x. The
error in norm(a*x-b) then should be within the machine precision, even
for ill-conditioned matrices.

Also, the SVD based least squares code fails completely for
non-square matrices in Julia 0.3.

Note that the pseudoinverse is not precise in Matlab/Octave either.
The above test should be accurate to a square root of machine precision
in this case.

The tests scripts and results for Julia 0.3, Octave 3.4.3, and Matlab R2012a follow:

---

Julia 0.3

After fixes:

julia> include("test_rdiv_lsquares.jl")
6.602160352047045e-14

julia> include("test_svd_lsquares.jl")
6.533966983681165e-14

julia> include("test_pinv.jl")
9.439277699982046e-7

---

Julia 0.3

Before fixes:

julia> include("test_rdiv_lsquares.jl")
8.788416928121552e-7

julia> include("test_svd_lsquares.jl")
ERROR: DimensionMismatch("assigned 100 elements to length 13 destination")
in assign_bool_vector_1d! at array.jl:382
in \ at linalg/factorization.jl:573
in include at ./boot.jl:242
in include_from_node1 at ./loading.jl:128
while loading /local/tmp/zng2/github-contrib/julia/test_svd_lsquares.jl, in expression starting on line 25

julia> include("test_pinv.jl")
0.00035320238415246054

---

Octave 3.4.3

octave> test_rdiv_lsquares
ans = 3.8725e-14

octave> test_svd_lsquares
ans = 3.9309e-14

octave> test_pinv
ans = 4.3197e-04

---

Matlab R2012a

matlab> test_rdiv_lsquares
ans =

1.0821e-13

matlab> test_svd_lsquares
ans =

1.0564e-13

octave> test_pinv
ans =

2.0777e-05

##########################################################################

## === test_rdiv_lsquares.jl ===

function hilb(n::Integer)
a=Array(typeof(1.0),n,n)
for i=1:n
for j=1:n
a[j,i]=1.0/(i+j-1.0)
end
end
return a
end

m = 1000
n = 100

# form an ill-conditioned matrix
a = randn(m,n) * hilb(n);

# create a right hand side in the range of a
x0 = randn(n);
b = a*x0;

# find the least squares solution
x = a\b;

# the error in rhs should be within machine precision
norm(a*x - b)

##########################################################################

%% === test_rdiv_lsquares.m ===

m = 1000
n = 100

% form an ill-conditioned matrix
a = randn(m,n) * hilb(n);

% create a right hand side in the range of a
x0 = randn(n,1);
b = a*x0;

% find the least squares solution
x = a\b;

% the error in rhs should be within machine precision
norm(a*x - b)

##########################################################################

## === test_svd_lsquares.jl ===

function hilb(n::Integer)
a=Array(typeof(1.0),n,n)
for i=1:n
for j=1:n
a[j,i]=1.0/(i+j-1.0)
end
end
return a
end

m = 1000
n = 100

# form an ill-conditioned matrix
a = randn(m,n) * hilb(n);

# create a right hand side in the range of a
x0 = randn(n);
b = a*x0;

# find the least squares solution
x = svdfact(a)\b;

# the error in rhs should be within machine precision
norm(a*x - b)

##########################################################################

%% === test_svd_lsquares.m ===

m = 1000
n = 100

% form an ill-conditioned matrix
a = randn(m,n) * hilb(n);

% create a right hand side in the range of a
x0 = randn(n,1);
b = a*x0;

% find the least squares solution
[u,s,v] = svd(a,0);
s = diag(s);
k = sum(s > eps*s(1))
s = diag(s);

s = s(1:k,1:k);
u = u(:,1:k);
v = v(:,1:k);
x = v*(s^(-1)*(u'*b));

% the error in rhs should be within machine precision
norm(a*x - b)

##########################################################################

## === test_pinv.jl ===

function hilb(n::Integer)
a=Array(typeof(1.0),n,n)
for i=1:n
for j=1:n
a[j,i]=1.0/(i+j-1.0)
end
end
return a
end

m = 1000
n = 100

# form an ill-conditioned matrix
a = randn(m,n) * hilb(n);

# create a right hand side in the range of a
x0 = randn(n);
b = a*x0;

# find the pseudo-inverse solution
x = pinv(a)*b;

# the error in rhs should be within the square root of machine precision
norm(a*x - b)

##########################################################################

%% === test_pinv.m ===

m = 1000
n = 100

% form an ill-conditioned matrix
a = randn(m,n) * hilb(n);

% create a right hand side in the range of a
x0 = randn(n,1);
b = a*x0;

% find the pseudo-inverse solution
x = pinv(a)*b;

% the error in rhs should be within the square root of machine precision
norm(a*x - b)

##########################################################################