Jul 012015

One of my goals when I started my website just after I finished grad school was to start sharing little snippets of code I had used or developed over the course of my undergraduate and graduate degree. One thing I did use, and still use frequently, is the Fourier transform. The Fourier transform is used in countlessly many algorithms in seismic imaging because of it’s relationship to sampling theory, ability to  identify frequency content in signals, and provide solutions to differential equations.

One day I stumbled upon a neat series of lectures by my undergraduate NSERC supervisor, Dr. Michael Lamoureux, for a Seismic Imaging Summer School in 2006. Hidden in a footnote of this set of lecture notes contains one of my favorite quotes about a mathematical theorem, namely the divergence theorem:

And, as I like to explain to my kids, the Divergence Theorem simply says you can measure how much methane gas a cow produces by either measuring how much is produced in each cubic centimeter of the cow, or simply by measuring how much leaks out its mouth and rear end, and other places on the surface of the cow. – Michael Lamoureux

However, I digress! In these lecture notes Michael describes a variety of theoretical considerations when investigating partial differential equations, and in particular the famous wave equation with initial conditions:

u_{tt} - c^2 \nabla^2 u = 0
u(x,0) = f(x)
u_{t}(x,0) = g(x)

In the third lecture of these notes, the Fourier transform is used to find a solution of the equation. This is also given as an example in Elias Stein’s famous book on harmonic analysis (though I have it written here as Page 395, I’m not sure which volume. It’s likely volume 1). In this case, the Fourier transform of the solution can be written as the sum of sines and cosines:

\hat{u}(k,t) = \hat{f}(k)cos(2\pi c|k|t) + \hat{g}(k)\frac{sin(2\pi c|k|)}{2\pi c|k|}

In the future, I might include the derivation though for now this will be sufficient. We simply apply the inverse Fourier transform to obtain a solution. In 2009 or so, I wrote a little Matlab code to demonstrate this solution. There are some scaling factors to match that of the Fourier transform implementation of Matlab. Luckily I’m a bit of a hoarder when it comes to old school work so here it is:

To be honest, there are a lot of poor programming practices going on here. For example when I wrote this, prior to my introduction into software development, I probably gave myself a nice little pat on the back for putting in those comments. Those who employ clean coding practices see where I’m going with this. If you take a look through the code, every time there is a comment it could be naturally refactored into a function with a descriptive name. Moreover, who am I trying to kid with variable names like “W”! No one would quickly realize that I’m referring to a grid of frequencies. Though it is academic code in nature, it’s still embarrassing!

Moving forward, we can then use this function to plot the solution for a couple different initial conditions:

Here are two sample videos:

Notice that we see some funky stuff happening at the boundaries. The mathematical solution given is valid on the entire plane, so it should just disappear off our screen. The fact that there is reflections at the boundary is an artifact of the discrete nature of computational mathematics, and these are typically handled using things like absorbing boundary conditions or perfectly matched layers. The easiest, but most computational inefficient way of handling it, would be to make an even larger grid and only record for the amount of time it takes for the reflection to reach the area you are recording.

Jan 242015

This caused me a bit of grief but the solution is somewhat anticlimactic for how long it took me! When I tried to mex compile some C++11 code with the “out of the box” mex compiler I ran into the following errors of the following type:

  1. warning: rvalue references are a C++11 extension [-Wc++11-extensions]
  2. warning: deleted function definitions are a C++11 extension [-Wc++11-extensions]
  3. no type named ‘shared_ptr’ in namespace ‘std’

And all sorts of weird stuff even though it compiled fine with clang++ on the command line. Essentially it came down to make the following adjustments to the mexopts file:


by adding adding -std=c++11 -stdlib=libc++ to CXXFLAGS:

And then forcing mex to use that file with -f and specifically including the path of the mex.h header file with -I:

Let’s see it in action with Lambda’s, introduced in C++11. Suppose the following is contained in lambda_hello_world.cpp.

Then inside MATLAB we execute the following:

Some references:

  1. http://www.cprogramming.com/c++11/c++11-lambda-closures.html
  2. http://stackoverflow.com/questions/16939734/matlab-mex-clang-c11-thread-undefined-symbols-error
  3. http://www.shawnlankton.com/2008/03/getting-started-with-mex-a-short-tutorial/
Jan 072015

So I’ve decided to take the plunge and have begun making repositories for my projects (old and new). I can’t believe I used to code without version control, it’s insanity! During grad school I had a neat project from my course in sparse matrix computation on Google Pagerank that I’m sharing on Github now. In my quest for making cool project names, I’ve dubbed this project LRPR(Low Rank Page Rank) and it can be found here: https://github.com/bee-rock/LRPR. Though it contains only one type of low rank approximation at the moment, perhaps that will change in the future.

The project report can be found here: The inner-outer solver for Pagerank combined with a subsampling scheme

The project starts out with a short introduction to Pagerank, formulated as en eigenvalue problem:

 Ax = x .


We investigate what the effect of a low rank approximation for the transition matrix has on the power method and an inner-outer iteration for solving the Pagerank problem. The purpose of the low rank approximation is two fold: (1) to reduce memory requirements (2) to decrease computational time. We show that we see an improvement in storage requirements and a decrease in computational time if we discard the time it takes to perform the low rank approximation, however at the sacrifice of accuracy.