Tags: #octave 

Boosting the Octave `hilb` function

Today, I read an interesting article about Hilbert Matrices in N. Higham’s blog and stumbled over these simple elegant two lines:

function H = hilb_higham (n)
  j = 1:n;
  H = 1 ./ (j' + j - 1);

I checked with the implementation of Octave. The currently implemented approach is not as elegant as it can be and includes a for-loop:

function H = hilb_octave (n)
  H = zeros (n);
  tmp = 1:n;
  for i = 1:n
    H(i, :) = 1.0 ./ tmp;
    tmp += 1;

Fun fact: The “MATLAB implicit expansion feature” mentioned in Higham’s article was introduced in Matlab R2016b. This feature is called broadcasting and was introduced in Octave 3.6.0 back in 2012.

Comparing the implementations clearly shows, that Higham’s solution should be taken.

In the figure generated by the test code below, smaller values are better.

N = 1000;
time_octave_hilb = zeros (1, N);
time_higham_hilb = zeros (1, N);

for i = 1:N
  hilb_octave (i);
  time_octave_hilb(i) = toc;

for i = 1:N
  hilb_higham (i);
  time_higham_hilb(i) = toc;

semilogy (1:N, time_octave_hilb, 'b');
hold on;
semilogy (1:N, time_higham_hilb, 'r');
legend ({"Octave 5.2.0", "Higham"});
xlabel ("Hilbert Matrix Dimension");
ylabel ("Time in seconds");


Octave uses Higham’s approach since cset 627da618dcc4.

July 01, 2020 (Version 1)

Download the Jupyter Notebook.

(C) 2017 — 2020 Kai Torben Ohlhus. This work is licensed under CC BY 4.0. Page design adapted from minima and researcher. Get the sources on GitHub.