Dan Monroe

Explorations in Math and Engineering

[when-im-alone-i-divide-by-zero]

Recursive Least Squares Stacking

My favorite part about linear algebra is taking an iterative function and stacking it into a beautiful vectorized algorithm.  In Recursive Least Squares you must choose a forgetting factor that defines how sensitive the algorithm is to changes in signal statistics, and in practice this usually amounts to a line search on that parameter.  In non-standard form, the update equations are

\[\begin{align} \bf{r}_{xx}(n) =& \sum_{k=0}^n f^{n-k}\bf{x}(k)\left[\bf{x}(k)\right]^T\cr \bf{r}_{dx}(n) =& \sum_{k=0}^n f^{n-k}d(k)\bf{x}(k)\cr \bf{h}(n) =& \bf{r}_x(n)^{-1}\bf{r}_{dx}(n)\cr y(n) =& \left[\bf{h}(n)\right]^T\bf{x}(n)\cr \end{align}\] where \(f\) is the forgetting factor, \(\bf{x}(n)\) is a \(L \times 1\) vector of filter inputs, \(d(n)\) is the desired signal, \(\bf{r}_{dx}(n)\) is the \(L \times 1\) cross-correlation vector, \(\bf{r}_{xx}(n)\) is the \(L \times L\) autocorrelation matrix of \(\bf{x}(n)\), \(\bf{h}(n)\) is the \(L \times 1\) coefficient solution, and \(y(n)\) is the scalar output (all varying with time \(n\)).  To avoid divison by \(0\), \(\bf{r}_{xx}(n)\) is initialized as \(\delta\bf{I}\) where \(\delta\) is a small constant.

Let us consider the case where \(L=1\).  We can begin by stacking all of the values of \(y(n)\) into an \(N\times 1\) matrix \(\bf{Y}\), which is just a convolution between \(\bf{x}(n)\) and \(\bf{h}(n)\).  We can define \(\bf{X}\) as a convolution matrix of \(\bf{x}(n)\), which is just

\[\bf{X} = \text{diag}(\bf{x}(n))\] and we can let \(\bf{H}_1\) be an \(N\times 1\) vector containing all values of \(\bf{h}_1(n)\) (remember, we are only using one coefficient), giving us

\[\bf{Y} = \bf{X}\bf{H}_1.\]

\(\bf{H}_1\) can be factored into

\[\bf{H}_1 = \bf{D_1}\bf{F}\bf{P},\]

where \(P\) is the inverse of  \(\text{diag}(\bf{R}_1\bf{F})\) and

\[\begin{align} \bf{D}_1 =& \left[\begin{array}{cccccc} &d(0)x_1(0) & 0 & 0 & \ldots & \cr &d(1)x_1(1) & d(0)x_1(0) & 0 &\ldots&\cr &d(2)x_1(2) & d(1)x_1(1) & d(0)x_1(0) &\ldots& \cr &\vdots & & & & \end{array}\right]\cr \bf{R}_1 =& \left[\begin{array}{cccccc} &x_1(0)^2 & \delta & 0 & \ldots & \cr &x_1(1)^2 & x_1(0)^2 & \delta &\ldots&\cr &x_1(2)^2 & x_1(1)^2 & x_1(0)^2 &\ldots& \cr &\vdots & & & & \end{array}\right]\cr \bf{F} =& \left[\begin{array}{ccccccc} &1 & f & f^2 & \ldots & f^n & \end{array}\right]^T\end{align}\]

giving us

\[\bf{Y} = \bf{X}\bf{D}_1\bf{F}\bf{P}.\]

This formulation gives us a way to alter the forgetting factor without recalculating all of the statistics.  This method is definitely less memory-efficient, and in a fast language like C this might be less processor efficient too, but in vectorized languages like MATLAB this reformulation increases speed by quite a bit, depending on the problem size and the size of your parameter space.  Sometimes it makes more sense to iterate with a for-loop, other times the vectorized method is best.  Here are some results computing the basic adaptive least squares solution with 196 different values of \(f\).  The for-loop algorithm is calculated using the recursive definitions of the statistics (because for filter length 1 it's faster than the more complicated realizations of the filter).

\[\begin{array}{ccc} \text{N (samples)} & \text{For-Loop} & \text{Vectorized} \cr  1000 & 3.097021 & 0.450168 \cr 2000 & 6.134475 & 1.986143 \cr 3000 & 9.314577 & 4.969422 \cr 4000 & 12.453965 & 9.253738 \cr 5000 & 15.420053 & 15.147522. \cr \end{array}\]

The numerical differences were, at worst, on the order of \(10^{-7}\).  This method works well for filters of length 1, and it even generalizes to filters with equality constraints, but that comprises an extremely small subset of all the possible filter lengths.  The "numerator" (everything except \(\bf{P}\)) can be generalized for an arbitrary length filter, giving

\[\bf{Y}_{\text{num}} = (\bf{X}_1\bf{D}_1 + \bf{X}_2\bf{D}_2 + \ldots + \bf{X}_L\bf{D}_L)\bf{F}\]

where \(\bf{X}_i = \text{diag}(\bf{x}(n-i+1))\) and

\[\bf{D}_i = \left[\begin{array}{cccccc} &d(0)x(0-i+1) & 0 & 0 & \ldots & \cr &d(1)x(1-i+1) & d(0)x(0-i+1) & 0 &\ldots&\cr &d(2)x(2-i+1) & d(1)x(1-i+1) & d(0)x(0-i+1) &\ldots& \cr &\vdots & & & & \end{array}\right]\]

assuming all initial conditions of \(\bf{x}(n)\) (for \( n < 0 \) ) are \(0\).

The issue comes with inverting the autocorrelation matrix at each time \(n\): apart from any strange linear algebra that I am ignorant of (maybe tensors) I have no way of compactly representing this operation.   A possible solution to this is pre-whitening the inputs: that is essentially the role of the autocorrelation inverse, but instead of whitening for every new sample and changing slightly with every statistical nuance of the data you will have a one-shot whitening step over the average statistics of all of the input data.  This is now no longer an exact approach, but will get a "ballpark" estimate.

[when-im-alone-i-divide-by-zero]

Estimating the Envelope of an Ultrasound Pulse

I have been doing some work with nondestructive evaluation, and one of the things I have been experimenting with is parameter estimation for an ultrasonic pulse model.  A pulse might look something like this:

We can characterize this as a gaussian-damped sunusoid, where

\[\begin{aligned} d(t) =& A\cos\left(\theta + 2\pi f_c (t-\tau) + \beta (t-\tau)^2\right)e^{-\alpha(t-\tau)^2}\cr =& \cos\left(\theta + 2\pi f_c (t-\tau) + \beta (t-\tau)^2\right)e^{\ln(A) + \alpha t^2 -2\alpha\tau t + \alpha\tau^2}\end{aligned}.\]

The Gaussian can be estimated in a couple of ways.  I started with an exponential curve fit of all maxima of the absolute value of the waveform (because the minimum peaks follow the same shape as the maximums).  The curve fit can be formulated as a least-norm optimization problem,

\[\begin{aligned} & \underset{\bf{h}}{\text{minimize}} & & \|\ln(\bf{|d|}) - \bf{X}\bf{h}\| \end{aligned}\]

where \(\ln(|\bf{d}|)\) is a vector of the absolute natural log of the signal, \(\bf{X}\) is a 2nd order Vandermonde matrix of our sample indices, and \(\bf{h}\) is a vector of the coefficients we are estimating, assuming our envelope is of the form

\[f_{\text{env}}(t) = e^{h_0 + h_1t + h_2t^2}.\]

The coefficients are related to our original paramters in the following way:

\[\begin{aligned} h_2 =& \alpha \cr h_1 =& -2\alpha \tau \cr h_0 =& \ln(A) + \alpha\tau^2. \end{aligned}.\]

The least-norm approach, even using the L1 norm, doesn't do a fantastic job of fitting the exponential.  The resulting curve is underfitted and misses the true envelope of the signal.

This underestimation is a side effect of regression on log-scaled values: less effort is expended minimizing errors for data points that are further away from zero because larger values become more and more compressed, while the values below 1 become more sensitive to errors.  We can mitigate this problem by shrinking our data window so more emphasis is placed on the higher amplitude points near \(\tau\).

It's still not perfect because of the way we are estimating the envelope of the signal.  If you notice, this fitted curve does exactly what we are asking it to: it passes through almost all of the peak values in the pulse, but it does not appear to rest on top of the waveform itself.  We can do much better by tweaking the way we pose the problem: we will add a constraint that the resulting curve is on or above the data, which should give us a true envelope.  It will also simplify preprocessing, getting us around the need to find the peaks of the waveform.  The new constrained problem is

\[\begin{aligned} & \underset{\bf{h}}{\text{minimize}} & & \|\ln(\bf{|d|}) - \bf{X}\bf{h}\| \cr & \text{subject to} & & \left(\ln(\bf{|d|}) - \bf{X}\bf{h}\right) \leq \bf{0} \end{aligned}\]

The new estimate fits much better.

Here are the parameter estimates, to full precision, of both methods (with the original window size):

\[\begin{array}{cccc} & \text{Original} & \text{Unconstrained} & \text{Constrained} \cr A & 72 & 63.44982 & 72.06581 \cr \tau & 102 & 103.0854 & 101.973 \cr \alpha & -0.003 & -0.003100467 & -0.002988858. \end{array}\]

They are both close, except the unconstrained method is pretty far off on estimating gain, \(A\).  Let's try this on some real data.  This is a 5Mhz backsurface echo sampled at 100MHz.

Applying the fitting techniques does not yield great results.

The problem is that there is some dissimilarity between the proposed mathematical model and what the pulse is actually doing.  The real pulse is asymmetric and doesn't fit into a true Gaussian shape, and beyond the main signal are small ripples present that distort the estimate.

We can first trim the pulse so we aren't looking outside of the main signal.  We can also relax the constraint in such a way that our solution will ignore the outside points and focus more on the center of the signal.  The original constraint says the fitted signal needs to lie on or above the original waveform; in other words, the error we are minimizing needs to be nonpositive.  If we relax the bound on the error, or allow it to become positive as we move away from the center of the pulse, we can more accurately model the part of the waveform we care about and ignore the rest.  Let us reformulate our problem and instead solve the following:

\[\begin{aligned} & \underset{\bf{h}}{\text{minimize}} & & \|\ln(\bf{|d|}) - \bf{X}\bf{h}\| \cr & \text{subject to} & & \left(\ln(\bf{|d|}) - \bf{X}\bf{h}\right) \leq \bf{F} \end{aligned}\]

where \(\bf{F}\) is a vector of constraints that relaxes as we move away from the peak of the function.  For this problem I defined the relaxation as follows:

\[\bf{F} = c|\max(\bf{d}) - \bf{d}|^k\]

where \(c\) and \(k\) are tuneable constants.  I set \(c = 0.05\) and \(k = .5\) and narrowed the window to eliminate the noise and focus on the most central portion of the pulse, allowing only 7 peaks near the center of the signal.  Narrowing that far in gave very similar estimates for both the constrained an unconstrained methods.

Zooming out, we can see the final approximation of the real ultrasound pulse.

It's not perfect but it's close, maybe the best we could do without changing the model.  This could have been entirely done by hand, tuning parameters to achieve the same results, but this type of method could also easily get a first estimate to begin working with.  If nothing else, this has been a fun little exercise in optimization.

[when-im-alone-i-divide-by-zero]

Another way to look at linear least squares

In least squares the objective is to estimate a set of observations given some inputs and a target signal.  Let our inputs and targets be zero-mean matrices where \(\bf{y}\) is a row vector of observations, \(\hat{\bf{y}} = \bf{H}^T\bf{X}\) is our linear model, \(\bf{H}\) is a column vector of coefficients, and \(\bf{X}\) is a matrix of inputs with each row being a different input signal.  We find the optimum Least Squares coefficients by minimizing the squared error between the target and the model, which can be written

\[\begin{align} e^2 &= \left(\bf{y}-\bf{H}^T\bf{X}\right)\left(\bf{y}-\bf{H}^T\bf{X}\right)^T\cr &= \bf{y}\bf{y}^T - 2\bf{H}^T\bf{X}\bf{y}^T + \bf{H}^T\bf{X}\bf{X}^T\bf{H}. \end{align}\]

The solution is well known and is expressed by the normal equations

\[\begin{align}\bf{X}\bf{X}^T\bf{H} =& \bf{X}\bf{y}^T \cr \bf{H} =& \left(\bf{X}\bf{X}^T\right)^{-1}\bf{X}\bf{y}^T, \end{align}\]

Then, the optimal estimate is \[\hat{\bf{y}} = \bf{X}^T\left(\bf{X}\bf{X}^T\right)^{-1}\bf{X}\bf{y}^T.\] The terms \(\bf{X}\bf{X}^T\) and \(\bf{X}\bf{y}^T\) are the sample autocovariance of \(\bf{X}\) and cross covariance between \(\bf{y}\) and \(\bf{X}\). If you factor the matrix inverse into a product of square roots you get \[\hat{\bf{y}} = \left(\left(\bf{X}\bf{X}^T\right)^{-1/2}\bf{X}\right)^T \left(\left(\bf{X}\bf{X}^T\right)^{-1/2}\bf{X}\right)\bf{y}^T.\]

The least squares solution is doing a few things here in cascade: it is performing a whitening transformation on the input \(\bf{X}\), and it is also computing the cross-covariance of each whitened \(\bf{X}\) and the target signal \(\bf{y}\) to generate each respective linear coefficient.  This means each least squares coefficient is a measure of how each (decorrelated) input varies in proportion to the output.  The input can be expressed as

\[\bf{W} = \left(\bf{X}\bf{X}^T\right)^{-1/2}\bf{X},\] the least squares coefficients are \[\bf{H}_{w} = \bf{W}\bf{y}^T,\] and the optimum estimate is \[\hat{\bf{y}} = \bf{H}^T_w\bf{W}.\]

[when-im-alone-i-divide-by-zero]

Geometric Series Matrix Inverse

Let the matrix \(\bf{S}\) be defined as follows:

\[\bf{S} = \bf{I} + \bf{A} + \bf{A}^2 + ... + \bf{A}^N \]

where \(\bf{A}^n = \prod_1^n\bf{A}\).  If you subtract \(\bf{AS}\) from \(\bf{S}\) you get

\[\begin{align} \bf{S}-\bf{AS} &= \bf{I} - \bf{A}^{N+1} \cr \left(\bf{I}-\bf{A}\right)\bf{S} &= \bf{I} - \bf{A}^{N+1}\end{align}.\]

If \(\left(\bf{I}-\bf{A}\right)\) is invertible then

\[\bf{S} = \left(\bf{I}-\bf{A}\right)^{-1}\left(\bf{I} - \bf{A}^{N+1}\right).\]

If \(N \to \infty\) and

\[\lim_{N \to \infty} \bf{A}^N = 0,\]

where \(\bf{0}\) is the zero matrix, then

\[\bf{S} = \left(\bf{I} - \bf{A}\right)^{-1} = \bf{I} + \bf{A} + \bf{A}^2 + ...\]

Let \(\bf{B} = \bf{I}-\bf{A}\), or \(\bf{A} = \bf{I} - \bf{B}\).  Then

\[\bf{B}^{-1} = \bf{I} + \left(\bf{I}-\bf{B}\right) + \left(\bf{I}-\bf{B}\right)^2 + ...\]

If \(\left(\bf{I} - \bf{B}\right)^n\) approaches \(\bf{0}\) as \(n \to \infty\) then the series approximates the matrix inverse.  If for some \(k\) the term \(\left(\bf{I} - \bf{B}\right)^k\) is nilpotent the solution is exact when expanded to the \(k\)th term.

[when-im-alone-i-divide-by-zero]

59:46

Neat.  I want to find that Tom Cover talk.  These examples are along those lines.

[when-im-alone-i-divide-by-zero]

First Post

Thanks to Christian Perfect and MathJax I finally have a nice way of using latex-style math symbols in a web page.  Lookie here:

\[y(x) = e^{\alpha\sqrt{x}}\]

I guess this means I can blog about math and engineering now.

I've been wanting to get away from MATLAB for a while, and over the past few days I have been experimenting with new tools.  So far Anaconda and julia are my favorites, and I started using a tool called the Sage Notebook, which is a great free scratchpad/CAS.

I really like julia because it's similar to MATLAB but it lets you get away with less (numeric type matching is not automatic).  Here is the source code and the plot for the equation above:

N1 = -1000; N2 = 1000 N = N2-N1 + 1; a = -0.3; x = linspace(N1,N2,N)*(1+0im); #Cast to a complex type y = exp(a*sqrt(x)); figure(1) plot(x,y) xlabel(L"$x$"),ylabel(L"$y(x)$") title(L"$e^{\alpha\sqrt{x}}$"*", "*L"$\alpha = $"*string(a)) grid(b=1)