JuliaLang/julia — pull request #9008 by zgimbutas (2014-11-14 17:43:30+00:00)
Title: Optimize pinv for diagonal matrices

In general, the pseudoinverse is determined to only half of the machine
precision for general dense ill-conditioned matrices. For diagonal and
certain other matrices (e.g., tridiagonal matrices), the pseudoinverse
can be determined more accurately. We try to detect if the matrix is
diagonal immediately after the call to pinv and invert the diagonal elements
only without invoking ths svdfact function which is both faster and
more accurate.

Since the default threshold can not be determined in advance for arbitrary
matrices, provide a user-selectable parameter tol and provide a simple
scheme for automatic threshold selection.

pinv(A) sets the threshold tol to eps(one(T))*maximum(size(A)).

pinv(A,tol) sets the user-selectable threshold tol.

pinv(A,0) attempts to determine the threshold tol by forming two pseudoinverses
for tol=eps(one(T))*maximum(size(A)) and tol=sqrt(eps(one(T)) and estimating
the error in the first Moore-Penrose relation a*pinv(a)*a = a. This can be done with very little overhead (since svdvals is an expensive call) while using
about twice as much intermediate memory.

The following is the test script and results:

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

function onediag(m::Integer, n::Integer)
a=zeros(typeof(1.0),m,n)
for i=1:min(n,m)
a[i,i]=1.0/(i^5)
end
return a
end

function tridiag(m::Integer, n::Integer)
a=zeros(typeof(1.0),m,n)
for i=1:min(n,m)
a[i,i]=1.0/(i^5)
end
for i=1:min(n,m)-1
a[i+1,i]=1.2/(i^5)
a[1,i+1]=1.2/(i^5)
end
return a
end

function test_pinv(a)

println("=== julia/matlab pinv, default tol=eps(1.0)*max(size(a)) ===")
@time apinv = pinv(a);

println(vecnorm(a*apinv*a - a))
println(vecnorm(apinv*a*apinv - apinv))
println(vecnorm((a*apinv)' - a*apinv))
println(vecnorm((apinv*a)' - apinv*a))
x0 = randn(n); b = a*x0; x = apinv*b;
println(vecnorm(a*x-b))

println("=== julia pinv, tol=sqrt(eps(1.0)) ===")
@time apinv = pinv(a,sqrt(eps(1.0)));

println(vecnorm(a*apinv*a - a))
println(vecnorm(apinv*a*apinv - apinv))
println(vecnorm((a*apinv)' - a*apinv))
println(vecnorm((apinv*a)' - apinv*a))
x0 = randn(n); b = a*x0; x = apinv*b;
println(vecnorm(a*x-b))

println("=== julia pinv, adaptive threshold ===")
@time apinv = pinv(a,0);

println(vecnorm(a*apinv*a - a))
println(vecnorm(apinv*a*apinv - apinv))
println(vecnorm((a*apinv)' - a*apinv))
println(vecnorm((apinv*a)' - apinv*a))
x0 = randn(n); b = a*x0; x = apinv*b;
println(vecnorm(a*x-b))

end

m = 100
n = 1000

# form an ill-conditioned matrix
println("\n--- dense ill-conditioned matrix ---\n");
a = randn(m,n) * hilb(n);
test_pinv(a)

println("\n--- diagonal matrix ---\n");
a = onediag(m,n);
test_pinv(a)

println("\n--- tri-diagonal matrix ---\n");
a = tridiag(m,n);
test_pinv(a)

We are testing the Moore-Penrose pseudoinverse criteria, all output numbers should be as small as possible:

--- dense ill-conditioned matrix ---

=== julia/matlab pinv, default tol=eps(1.0)*max(size(a)) ===
elapsed time: 0.024110237 seconds (2979592 bytes allocated)
8.132311486154941e-5
202456.53269313014
8.717569398698316e-5
0.00042487526397591465
8.958400718585145e-6
=== julia pinv, tol=sqrt(eps(1.0)) ===
elapsed time: 0.026959505 seconds (2979464 bytes allocated)
1.1318678252271499e-7
0.00044092008873858384
8.220588998363644e-9
3.41201657124137e-8
1.9379690006220722e-8
=== julia pinv, adaptive threshold ===
elapsed time: 0.027885864 seconds (7306016 bytes allocated)
1.1318678252271499e-7
0.00044092008873858384
8.220588998363644e-9
3.41201657124137e-8
5.994505957502731e-8

--- diagonal matrix ---

=== julia/matlab pinv, default tol=eps(1.0)*max(size(a)) ===
elapsed time: 0.000258462 seconds (800048 bytes allocated)
5.423670138152395e-20
1.639026212099123e-6
0.0
0.0
1.694331095701549e-21
=== julia pinv, tol=sqrt(eps(1.0)) ===
elapsed time: 0.000260244 seconds (800048 bytes allocated)
5.423670138152395e-20
1.639026212099123e-6
0.0
0.0
5.43175421391522e-20
=== julia pinv, adaptive threshold ===
elapsed time: 0.000233371 seconds (800048 bytes allocated)
5.423670138152395e-20
1.639026212099123e-6
0.0
0.0
2.71084145559175e-20

--- tri-diagonal matrix ---

=== julia/matlab pinv, default tol=eps(1.0)*max(size(a)) ===
elapsed time: 0.009671382 seconds (2980904 bytes allocated)
2.6789350069963263e-15
5.6684013753107605
1.2426503022248955e-8
4.506044543963786e-5
2.0526552865923344e-15
=== julia pinv, tol=sqrt(eps(1.0)) ===
elapsed time: 0.009924386 seconds (2979752 bytes allocated)
5.000380818672045e-8
0.0019228562453838195
2.2514327087865425e-10
2.4024973745285115e-9
5.212374715982636e-8
=== julia pinv, adaptive threshold ===
elapsed time: 0.013035343 seconds (7307616 bytes allocated)
2.6789350069963263e-15
5.6684013753107605
1.2426503022248955e-8
4.506044543963786e-5
1.3740389626696697e-15