Vectorize Matrix Slice Extraction

There was an interesting thread on the Octave mailing lists recently about vectorizing matrix slice extraction. That is avoiding the dreaded for-loop in something like:

%% create a random sparse matrix
n = 10000; %% matrix dimension
nz = 300; %% number of non-zero elements

ia = floor(rand(nz,1)*n); %% dimension 1 indices
ja = floor(rand(nz,1)*n); %% dimension 2 indices
A(:,:) = sparse(ia, ja, rand(nz,1), n, n, n);

%% taking slices with a for loop, slow, poor scaling to large matrices:
for jj = nz:-1:1
slow_columns(:,jj) = A(:,ja(jj));
endfor

%% example using built-in sub2ind, fast, good scaling to large matrices:
jj = sub2ind(size(A)(1),ja);
fast_columns = A(:,jj);

This is a simplified example, the one in the thread shows how you can access a 5-dimensional matrix with 4 indexes. Knowing the built-in functions that come with your high-level language is well worth it, even if it means you have to do weird things to stay in the high-level language and still get reasonable speed. This allows you to avoid slow looping and use the quick compiled functions that someone else has already spent time and effort developing.

Unblurring Images

Tikhonov regularization is a pretty cool concept when applied to image processing. It lets you 'de-blur' an image, which is an ill-posed problem. The results won't be as spectacular as CSI, but still cool nonetheless.

Original

The original image, CSI fans pretend it is a license plate or a criminal's face in a low quality surveillance video, was made by taking a screenshot of some words in a word processor. The image doesn't have to be character based, but characters make good practice problems.

Blurring

We then blur the image with a Toeplitz matrix. In Octave that's pretty easy:

im = imread("foo-fiddle.png");
[m, n] = size(im);

sigma_true = 10;
p = exp(-(0:n-1).^2 / (2 * sigma_true^2))';
p = p + [0; p(n:-1:2)];
p = p / sum(p);

G_true = toeplitz(p);

im = G_true*double(im)*G_true';

Generally after doing some floating point operations on an image you won't end up with something still scaled between 0 and 255, so the little function shown below rescales a vector's values to the proper range for a grayscale image.

function im = image_rescale(double_vector)
%% scale values linearly to lie between 0 and 255
m = max(double_vector);
n = min(double_vector);
scale = max(double_vector) - min(double_vector);
im = double_vector - (min(double_vector) + scale/2);
im = im ./ scale;
im = (255 .* im) + 127.5;
endfunction

De-Blurring

This is where the regularization comes in. Since the problem is poorly conditioned we can't just recover the image by straightforward inversion. We'll use the singular value decomposition of our blurring matrix to create a regularized inverse, this is similar to the pseudoinverse. If the singular value decomposition of the blurring matrix is


Then the regularized inverse is


Where the entries in the diagonal matrix are given by


So as alpha approaches zero we approach the un-regularized problem. Chapter 4 of Watkins has good coverage of doing SVDs in Matlab. The Octave to perform this process is shown below in the function tik_reg().

function im_tik = tik_reg(im,sigma,alpha)
[m,n] = size(im);
p = exp(-(0:n-1).^2 / (2 * sigma^2))';
p = p + [0; p(n:-1:2)];
p = p / sum(p);

G = toeplitz(p);

[U, S, V] = svd(G,0);

Sinv = S / (S^2 + alpha*eye(length(S),length(S))^2);

im_tik = V*Sinv*U'*double(im)*V*Sinv*U';
endfunction

This gives us two parameters (sigma and alpha) we have to guess in a real image processing problem (because we won't be the ones blurring the original). Another weakness of the method is that we assumed we knew the model that produced the blur, if the real physical process that caused the blur is not similar to Gaussian noise, then we won't be able to extract the original image.

Here's a couple for-loops in Octave that try several combinations of alpha and sigma.

sigma = [9,10,11];
alpha = [1e-6,1e-7,1e-8];

for j = 1:length(sigma)
for i = 1:length(alpha)
unblur = tik_reg(im, sigma(j), alpha(i));
unblur = reshape(image_rescale(reshape(unblur,m*n,1)),m,n);
fname = sprintf("unblur_%d%d.png",i,j);
imwrite(fname, uint8(unblur));
endfor
endfor

Optical character recognition probably wouldn't work on the output, but a person can probably decipher what the original image said.

CFD Integrity

Quantifying uncertainty in simulation predictions is something at which most modelling communities fall short, and the Computational Fluid Dynamics community is no exception. There is plenty published on uncertainty estimation, but the state of the practice lags the state of the art considerably in this area.

It does not help matters that it is very hard and expensive to do a thorough job of quantifying the uncertainty of the result of a complex calculation. If it takes a week to get one data point, doing a Monte Carlo sensitivity analysis is probably out of the question, and differentiating your code could take years. It also doesn't help that many of the decision makers who are the consumers of the CFD product don't really want to hear how uncertain the results are, often they want pretty pictures and a warm fuzzy feeling (colorful fluid dynamics is great for marketing too).

I highly recommend the article on quantifying uncertainty by P.J. Roache. It is well written, and it is a special treat to read a CFD paper with a sentence like,

Clearly, we are not interested in such worthless semantics and effete philosophizing, but in practical definitions, applied in the context of engineering and science accuracy.

Roache is addressing the hand wringing made popular in the "soft" academic communities about the ability to define "truth".

Other than being a good basic engineering practice, making sure users of your results are well aware of the uncertainty is a matter of integrity too. There's an obligation that comes with expertise. Boole's words on the subject still ring true. More recently Feynman said much the same thing about integrity as an expert advising others:

I'm talking about a specific, extra type of integrity that is not lying, but bending over backwards to show how you are maybe wrong, that you ought to have when acting as a scientist. And this is our responsibility as scientists, certainly to other scientists, and I think to laymen.

Convective Flux Jacobians

The Euler equations are an important system of equations in computational fluid dynamics (CFD). They are basically the full continuum equations of motion without the viscosity. Adding in the viscosity gives you the more realistic Navier-Stokes equations.

The conserved variable vector is

And the flux vector is

The flux Jacobian is simply the partial derivative of the flux with respect to the conserved variables (mass, momentum, energy).

The maxima to calculate the Jacobian is straightforward:

U : matrix([u1], [u2], [u3]);

f1 : u2;
f2 : u2^2 / u1 + p_exp;
f3 : (u2 / u1) * (u3 + p_exp);
F : matrix([f1], [f2], [f3]);

A : zeromatrix(nd, nd);

for j : 1 thru nd do
(
for i : 1 thru nd do
(
A[i,j] : diff(F[i,1], U[j,1])
)
);

Diagonalization

Using the chain rule lets you write the Euler equations in a quasi-linear form in terms of the flux Jacobian, now called A.

The flux Jacobian can be diagonalized into a matrix of right eigenvectors, diagonal matrix of eigenvalues, and a matrix of left eigenvectors.

The code to accomplish this in maxima is a little more involved, because of some substitutions needed to get everything in terms of u,gamma,a.

load("eigen");
VLV : eigenvectors(simp_A);

a_exp : (g * (g - 1) / u1) * (u3 - u2^2 / (2 * u1));
a_exp : sqrt(a_exp);
simp_a : ratsimp(ev(a_exp, u1 = rho, u2 = rho * u, u3 = rho * E));
a_eqn : a^2 = simp_a^2;
c_exp : g*(2*E-u^2);
a_eqn : facsum(a_eqn, c_exp);

/* put the eigenvectors into a matrix: */
V : columnvector( VLV[2] );
V : addcol( V,VLV[3],VLV[4] );
V : ratsimp( V );

/* put the eigenvalues into a diagonal matrix: */
load("diag");
L : diag( VLV[1][1] );

radsubstflag : true;
/* substitute in the speed of sound: */
V : ratsimp(ratsubst(a, simp_a, ratsimp(V)));
V : factor(ratsimp(ratsubst(a, simp_a, rootscontract(V))));
/* factor out gammas */
V : facsum(V,g);
c_a : part(solve(ratsubst(c,c_exp,a_eqn),c)[1],2);
V : ratsubst(c,c_exp,V);
V : ratsubst(ev(c_a),c,V);
V : facsum(V,u);

/* substitute in the speed of sound: */
L : ratsimp(ratsubst(a, simp_a, ratsimp(L)));

/* left eigenvectors: */
Vinv : factor( ratsimp( invert( V ) ) );
/* factor out gammas */
Vinv : facsum(Vinv,g);
d_exp : (g-1)/a^2;
Vinv : ratsubst(d,d_exp,Vinv);

Conclusion

The thing to note about this result is that there are three distinct wave speeds in the one dimensional Euler equations, u-a, u+a, and u, where a is the speed of sound. There are lots of possible ways to write the eigenvectors, in terms of the enthalpy, or the energy, or the pressure, but the form above with only the flow velocity, the speed of sound, and the ratio of specific heats is pretty elegant. It really illustrates what a fundamental quantity the speed of sound is for this set of equations.

Chapter 3 of Riemann Solvers and Numerical Methods for Fluid Dynamics: A Practical Introduction has a good work-up of this development and shows some alternative formulations for the flux Jacobian.

Kronecker Product of Sparse Matrices

Gilbert Strang discusses the use of a neat built-in function in Octave, kron(), that's pretty useful for doing multidimensional finite difference methods on PDEs:


kron() takes the Kroenecker product of two matrices. This is especially useful when constructing the large, sparse matrices needed for finite difference approximations.

K is the tridiagonal matrix created for a central difference approximation to the second derivative:

n = 5;
K = spdiag(2*ones(n,1),0) + spdiag(-ones(n-1,1),-1) + spdiag(-ones(n-1,1),1);
spy(K)

K is the differences for a 1D row of points, now we can create the 2-dimensional version for a plane of points with kron():

K2D = kron(speye(n),K) + kron(K,speye(n));
spy(K2D);

The extension to 3 dimensions is exactly as before, now we have a matrix describing the differences between points on the same line in the first sub and superdiagonal, differences across lines in the next bands out from those, and differences across planes in the furthest out bands:

K3D = kron(speye(n),K2D) + kron(K2D,speye(n));
spy(K3D);

Pretty simple! Now you're ready to solve multidimensional partial differential equations.

Fibonacci Numbers

Gilbert Strang has an interesting example in one of his linear algebra lectures on eigenvalues that is relevant to the articles I've written on speeding up calculations in Octave.

You might be familiar with the Fibonacci numbers, they are calculated by a two term recurrence:


We've already seen that recasting a calculation in terms of a matrix vector multiply or matrix vector solve is an important skill for writing good fast code in Octave. Strang demonstrates this at about halfway through the above lecture. The first order system that describes the generation of Fibonacci numbers is:


The nth term of the sequence is calculated by taking the nth power of A and applying it to the initial conditions:


When you diagonalize a matrix the powers of the matrix become very easy to calculate:

If you aren't familiar with the concept of eigenvalues and eigenvectors and diagonalizing a matrix I would highly recommend the linear algebra lectures on MIT's OCW website.

Here's the Octave code to accomplish the above. If you want, you can use tic() and toc() for timing. On my laptop the matrix method is about three times faster than the for loop method. What kind of results do you get? This advantage will grow as you go further into the sequence because the amount of work stays constant for the matrix method, but grows linearly for the recurrence method.

%% Our 'fibonacci' operator:
A = [1,1;1,0]
%% our initial conditions:
u_0 = [0;1];

%% get the matrix of eignevectors, V, and the matrix of eigenvalues, lambda
[V,lambda] = eig(A)

%% A = V lambda V^{-1}

Vinv = inverse(V)

%% This gives us an idea of the error introduced by calculating the
%% inverse with finite precision arithmetic:
Error = A - V*lambda*Vinv

n = 1400;
F = zeros(n,2);
%% Find the 100th Fibonaci number:
F(100,:) = V * (lambda^100) * Vinv * u_0;

%% check our answer with a simple loop:
fibo = zeros(n,1);
fibo(1) = 0;
fibo(2) = 1;
for i=3:n
fibo(i) = fibo(i-1) + fibo(i-2);
endfor

for i=2:n
F(i,:) = V * (lambda^i) * Vinv *u_0;
endfor

relative_error = abs( ( F(3:n,2) - fibo(3:n) ) ./ fibo(3:n) );

One thing to note is that one of the eigenvalues of our matrix is greater than one. This has implications for the numerical stability of the above method for calculating terms in the Fibonacci sequence. It means that the initial tiny round-off error in our eigenvector matrix inverse will be amplified as we continue to calculate higher and higher terms in the sequence. This behaviour is shown in the plot below, luckily the error starts out very small, and doesn't grow too large before we run out of range in the 32bit floating point numbers on the machine.

What's the answer? The 100th Fibonacci number is approximately 2.1892e20.

How to Speed Up Octave

Speed seems like a continuing theme on the Octave mailing lists, a recent thread on the list of "comparing the executive speed with Matlab" demonstrates that there really is very little speed overhead when you use Octave smartly. If you read down far enough, it also demonstrates that the Intel Fortran compiler is awesome!

There are four basic methods of speeding up a chunk of Octave code (in order of increasing speed):

1. Allocate memory smartly (using for loops)


2. Vectorize your code (avoiding for loops)
   


3. Use sparse matrix operators in place of for loops you can't vectorize


4. Write a compiled function in C/C++/Fortran to call from Octave

The first three options keep you in pure Octave. The fourth option, while requiring you to use a compiled language is really not too bad because you can use all of the types (vectors, matrices, etc) defined in the Octave libraries to write your function. That's almost as easy as doing it in Octave.

Speeding Up For Loops

This one really amounts to not being foolish. The only way to speed up for loops is to make sure you allocate all of your memory (using zeros() for instance) before you enter the for loop. The Octave Manual covers this topic here. The better thing to do is vectorize your code and avoid for loops all together.

Vectorizing

What you really want to do is cast your calculation as a series of vector and matrix operations, or use the built-in (usually compiled) functions that operate directly on vectors. The post Is Octave Slow? demonstrates the sorts of speed up you can get when you get rid of for loops.

Sparse Matrices for Data Dependancy

What about for loops you can't vectorize? The ones where the result at iteration n depends on what you calculated in iteration n-1, n-2, etc. You can define a sparse matrix that represents your for loop and then do a sparse matrix solve. This is much faster because you are now using compiled functions. The drawback is that you are trading some increase (sometimes substantial) in memory usage to buy this speed up.

• Is Octave Slow? covers this approach as well.


• More on Speeding up Octave covers the details on sparse matrix allocation.

Compile It

The manual has a good intro on using mkoctfile to compile your function and make it callable from Octave. Probably the most useful thing is that you can use the Octave types, so you don't have to code up your own matrix or vector type. At this point the speed really has little to do with Octave and everything to do with the quality of the compiler you use and how efficiently you code up your calculation.

Conclusion

The real power of Octave isn't that you can get your calculations to be nearly as fast as if you wrote them all in a compiled language. It's that you can prototype a calculation quickly and make sure that it is correct, and then there is a fairly well defined path to follow to speed up your correct calculation.

Thresholding for Edge Detection

The real world is noisy, and the differencing operations of edge detection tend to magnify that noise. We can apply some statistical methods to our edge detection to hopefully overcome that. A very basic approach is to threshold based on a histogram of the gradient magnitude. Taking a look at the distribution of gradient magnitude for our trusty old electricity meter (shown below), it looks like a pixel level below about 235 puts us safely out into the tail of the distribution. Most pixels (in this image) are not edge pixels, so they are clustered on the high end of the scale.

This is pretty easy to do with numpy and Image (PIL):

import numpy
import Image
pix = numpy.array( img.getdata() )
for i in xrange(0,pix.size):
if pix[i] < 235: pix[i] = 0 else: pix[i] = 255 img_out.putdata( pix )

Of course the histogram shape will depend on how much noise and how many edges are in the particular picture. This one had black dials on a white background, and not much noise so thresholding is a pretty fruitful strategy for really pulling out the edges.

You can make an intuitive argument that gradient magnitude distributions will have roughly similar shapes to that shown above. Especially as image resolutions steadily increase, the vast majority of pixels will not be on an edge. Sort of the curse of dimensionality in reverse, I guess that makes it a blessing.

Incompatibility between Matplotlib and Numpy

If you are running Fedora 9 you may have run in to this problem (at least one other guy did). Taking a look at the bugzilla over at RedHat shows that this has been a recurring problem over the past couple years (first cropping up in FC5). It seems that the packagers for numpy and matplotlib are having trouble synchronizing the updates for numpy and matplotlib. They are fast moving projects so one occasionally changes the API quietly and breaks the other. Of course there are other dependencies that probably play in to the decision of when to update a particular package, but this is the one that breaks my workflow.

The easiest fix for me was to download the latest matplotlib sources and compile/install them manually:

$ tar -zxvf matplotlib-0.98.5.2.tar.gz
$ cd matplotlib-0.98.5.2
$ python setup.py build
$ su
$ python setup.py install

This is less than ideal, because now I'm not running a "stock" Fedora where all of my software updates are taken care of automagically, but the increase in work-load is negligible.

A work-around for the incompatibility between numpy and matplotlib that requires editing some of the python sources:
http://www.mail-archive.com/numpy-discussion@scipy.org/msg13753.html

Or you can just upgrade to Fedora 10. Unfortunately that's not an option for me because of a closed proprietary device driver that relies on < 2.6.27 kernels. I haven't compiled my own kernels since FC1 days, and don't really feel like getting back into the habit.