Simple linear inversion [code]

A couple of days ago I was asked to make a figure showing a simple linear mixing/inversion problem. It's a very straightforward task, but I thought this would be a good opportunity to share some MATLAB code online for the first time.

You can find the code here. First I'll explain what it produces, and then I'll talk a little about the code itself.

My task was to show the basics of spectral mixing. In spectroscopy, light is used to distinguish different chemicals by which wavelengths they absorb or reflect. A spectrum looks something like this:

The horizontal axis is in wavenumbers, which is just the inverse of the wavelength, $k = 1/\lambda$. The vertical axis gives the percentage of light absorbed if shone through the test sample. If there's a peak, it means that the molecules in the sample can be excited by that wavelength, which is a clue as to what kind of chemical it is.

Here I'm just making up spectra: there are no chemicals that correspond to these Gaussian lineshapes, but I needed something I could do quickly, and Gaussians are quick. I wanted to create a sample with several chemical contributions drawn from a large library of chemicals. Each member of the library is assigned a Gaussian spectrum with a randomly generated height, center, and width. The library is called $\mathbf{A}$ and looks like this:

I randomly select a small number of these, say five. Then I have to decide how much of each chemical is in the simulated sample. I do this by choosing a coefficient for each spectrum randomly, and then normalizing to make sure the coefficients add up to 1. That is, I choose a percent contribution for each spectrum in the sample. If we represent this so-called density vector is by $D$, then the composite spectrum is just the product $S=\mathbf{A}D$. It looks like this:

The black line shows the composite spectru, and the colored lines show the individual spectra multiplied by their percent contributions.

So far, we've looked at the forward model. That is, you tell me what the sample is made of, and I can tell you what the measurement should look like. This model is really simple. It just says that you take a weighted sum of the individual spectra. It's a linear system, basically meaning that if you double the inputs you get double the output.

Now we want to do the inverse problem: you tell me what the measurement looks like, and I'll tell you what the sample is made of. In general, inverse problems are close to impossible. Imagine if I told you that I weighed an object and discovered that it weighed one pound, and I ask you to identify the object. You might be able to make a good guess, but you don't have enough information to solve the problem exactly. I'll talk a lot more about making decisions under uncertainty, but the point for now is that you can't expect all inverse problems to have unique answers.

This problem is different. Not only do we have enough information to solve it, but inverting linear problems can be done with one line of simple code. The density is recovered by the equation
\begin{equation}
D = \mathbf{A}^{-1}S
\end{equation}
where  $\mathbf{A}^{-1}$ is the inverse of $\mathbf{A}$. We could find it by hand, but MATLAB does it for us if we call the function pinv().

The result is a density for each member of the spectral library, most of which are zero since they're not in the sample. In fact, we can quickly check which members are non-zero (actually about a threshold, just in case of machine errors). We can compare the recovered members to the members we chose to make sure the inversion worked.

***
Notes on the code

The main code is spectral_mixing, and there are three supporting functions, create-library, plot_spectsum, and plot_composite_spectrum. These output nicer plots that MATLAB does natively.

If you want to export these graphics as .eps files, I recommend using the print2eps function included in the wonderful export_fig package by Yair Altman.