Ganesh Swami

A couple of months back, Prof. Gilbert Strang had come to SFU as part of the distinguished speakers series (as I had written earlier.) After the talk, I got a chance to chat with the guru. I owe a lot of my understanding of linear algebra to his books and his kick-ass animations (check out these eigen-analysis demos,) and so I asked him about the single most important topic in linear algebra. Without hesitation, he immediately responded: “the SVD!”

The Singular Value Decomposition is the swiss-knife of linear algebra. Every matrix Y can be factored into three matrices: U, S, and V as

U and V are orthogonal matrices and S is a diagonal matrix. Some uses of this factorization of the matrix Y: calculating 2-norms, Frobenius norms, ranks, null spaces and ranges, (pseudo)inverses and determinants, and by extension solving systems of equations (exact, over- and under- determined), eigenvalues and eigenvectors, and approximations.

Almost every problem in engineering becomes an optimization problem (as in reducing the error under some norm) and the method of least squares makes extensive use of the SVD.

The stage has been set for the next few posts of mine.

I’ve recently become a heavy user of the numerical capabilities of Python. I’ve written about my experiments before, but now I’m writing production quality code with numpy and matplotlib.

The above is an actual plot that I created for some Hall measurements I was doing. I was supposed for find functional relationships between temperature and majority charge carriers, which in my case were electrons because of the n-type doping. The simple case was a least squares fit: scipy.optimize.leastsq to the rescue. The more complicated part was solving a non-linear equation for roots and then doing a least squares fit. The root-finding module in scientific python provides lots of options. At this point, I can confidently say that this environment has more features than Octave.

Just today, I wanted to use the Fourier method on a differential equation (plug: the advantages of which are here) and numerical python with fft, fftshift and fftfreq are exact substitutes for their Matlab equivalents. You can also put actual LaTeX equations on plots, which is a major plus.

That is all.

I strongly think that unit tests are absolutely necessary to keep the overall quality of ones code high. I’ve written about Unit Tests before, so I’ll continue on that trend.

Most high performance numerical codes use some variant of Fortran or just raw C (which is my staple language.) In that case, as Micheal Feathers rightly points out, unit tests can “collide” with C.

I’ll want to expand on his solution of link seams.

GNU/Linux (and other Unix-like operating systems I think), provide a mechanism of loading shared object files at runtime. Using dlopen(3) and dlsym(3) you can re-create much of the reflection functionality that the Java community enjoys. So by loading your unit test dynamically at runtime, and having a predetermined symbol that behaves as the entry point, you can nicely wrap your module/function with a unit test. Interface this with ctest from the cmake project, and you are set. If you want to have a private function whose symbol isn’t visible on the outside, make it static. The nm utility will give you a list of public symbols in an object.

I’m pretty sure this is supported across the board. I remember doing this under DOS using DJGPP (remember that?) so this isn’t exactly novel.

For some work I was supposed to do later this term, I wanted a few synthetic vector/tensor fields. Something more complicated than a regular planar field.

I’ve used radial and tangential fields in electromagnetics (think solenoids and inductors) countless times, so it should be trivial to generate one, shouldn’t it? Unfortunately, I was getting mixed up in the minus signs somewhere and for the life of me couldn’t find out where. I checked and rechecked my math. I checked the usage of atan2(3) instead of atan(3). Finally, I had to resort to generating them by rotating the Cartesian basis and regenerating the tensors from the spectral components. Hacky, cludgy and doesn’t follow the DRY principle but this is just a test case. Yay tangential field:

To check the robustness of my algorithms later on, I need to add some noise to the fields and see how well they perform. Plagued by partial volume effects, diffusion tensor data are inherently very noisy, so it’ll be good to include noise as part of the algorithm development process. Right now, modeling noise in tensors is a very complicated process because the tensors themselves are built through a linear regression from diffusion weighted images. The noise is definitely not gaussian. Gaussian noise is easy, and that’s what I’m doing now until I fully understand how noise is transformed through the regression.

Another matter of complication is that regularizing (or denoising/smoothing) these fields is an active area of research. Extensions to the standard anisotropic edge-preserving filters like Perona-Malik are non-trivial (at least to me.) One of the technicalities of diffusion tensors is that they are positive definite. The positive definiteness is a physical manifestation as diffusion can only be zero at absolute zero (a great story for sci-fi.) I picked one of the many extensions just to test my workflow and here are the results. Pretty good I’d say for all the complications.

Lastly, opensource visualization tools seem to trip over the simplest of tasks. MayaVi seems to be the only one that can update the VTK pipeline when the source data changes. In other programs I have to rebuild the pipeline each time. As you can imagine, this gets tiring really fast.

I must be on a roll here. For the last few days, I’ve been spitting out about 500 words of coherent, easy to understand text for my coop report. I’ve never enjoyed writing technical documents so much before.

Everytime I intend to write a report, I always start with the LaTeX template that I made. This time, I decided to push the limits of what LaTeX has to offer.

I’ve included some pictures from my half-baked document. Includes pictures of Table of Contents, PDF indexes, Table of Notation, Glossary, and References. All of them cross-referenced and hyper-linked automatically. Click on the images for higher resolution. Resolution is not meant to be able to read the text, it’s still a bunch of FIXMEs.