Re: [Gretl-devel] basic Julia support
On Tue, 19 Jan 2016, Riccardo (Jack) Lucchetti wrote:
On Sat, 16 Jan 2016, Allin Cottrell wrote:
> My test was admittedly kinda silly, in that there's not really
> any reason to delegate to a "foreign" program stuff that gretl
> handles well natively. One would be more likely to get Julia to
> do MCMC or the like.
[snip]
This is probably more akin to a real-life case:
[case: delegating matrix inversion to julia]
I'm not sure we have a "real-life" case yet. I _am_ sure that there
exist problems that gretl could usefully send to julia for
substantial time-saving, but I don't think that matrix inversion is
one such. Consider the following, which mutiplies Jack's example up
to the point where the timings can be distinguished:
<hansl>
set echo off
set messages off
dim = 12
matrix X = mnormal(dim, dim)
set stopwatch
loop i=1..10000 -q
iX = inv(X)
endloop
printf "native inversion: %.3f secs\n\n", $stopwatch
mwrite(X, "X.mat")
set stopwatch
foreign language=julia
X = gretl_loadmat("X.mat", 0)
function f(X)
for i=1:10000
jiX = inv(X)
end
end
@time f(X)
gretl_export(inv(X), "invX.mat", 0)
end foreign
printf "delegated inversion: %.3f secs (with I/0 overhead)\n\n",
$stopwatch
check = X * iX - I(dim)
printf "gretl: max error = %g\n", max(abs(vec(check)))
jiX = mread("invX.mat")
check = X * jiX - I(dim)
printf "julia: max error = %g\n", max(abs(vec(check)))
</hansl>
On my intel sandybridge desktop I get the following:
<output>
native inversion: 0.323 secs
0.421968 seconds (313.88 k allocations: 96.529 MB, 7.02% gc time)
delegated inversion: 2.201 secs (with I/0 overhead)
gretl: max error = 1.98469e-15
julia: max error = 2.15561e-15
</output>
The "delegated" inversion time using julia is about 7 times greater
than the the native time. This comparison is obviously somewhat
unfair to julia since it includes set-up and I/O time for julia but
not gretl. Nonetheless, the internal julia time (indented a couple
of spaces) is close to the native time (in this run, a little
greater), so why delegate?
We're both using openblas for this sort of operation, so I see
little or no advantage in calling a "third-party" program. By the
way, the "max error" figures (relative to an identity matrix when
calculating A * inv(A)) are in both cases basically random numbers
on the verge of double-precision zero; gretl happens to come out
slightly better in the output above, but that's arbitrary.
Perhaps we should offer a prize for the first example of actual
time-saving (and/or improvement in accuracy) in sending a
calculation of the sort that might be required in gretl out to julia
for computation.
Allin