Fitting data to sum of exponentials

1

I would like to fit experimental data with exponential decay to an expression that looks like the following: f(t) = a_1*exp(-t/\tau_1)+a_2exp(-t/\tau_2) + ...+ a_n*exp(-t/\tau_n) , where the values of \tau_n are predetermined. Therefore, I would like to do a least square fit. I have a 1000 data points with 1000 corresponding values of t.

Given this, I had the following questions:

  1. If I want to do a least square fit with the form minimize |Ax-b|_2^2 , how would I structure A given that A is a function of t. If there were only one t value and a 100 values of \tau , then I could see that A would be a 100X100 matrix. What would it be in the case where I am fitting over a 1000 points of t? If I formulated this as a 3 dimentional matrix, it would be a 100X100X1000 matrix. However, I don’t think cvx supports 3D matrices.

  2. Would I use the sum_square(A*x-b) function in cvx?

Thank you in advance for your help!

2

Your model as presented above, is just unconstrained ordinary linear least squares (it is linear because your least squares (regression) equation is linear in the unknown coefficients). I am not providing judgment as to the merits, or lack thereof, of your formulation as ordinary linear least squares.

Let m be the number of data points. You can minimize

norm(A*x-b)

where

x is the column vector [a1;a2; ... ;an]  to be optimized

and

A is m by n, with one row per data point ti, consisting of
[exp(-ti/tau1) exp(-ti/tau2) ... exp(-ti/taun)]

and

  b is the m by 1 column vector of values of f(x) (i.e., the Left Hand Side).

In your example above, 100 values of tau corresponds to n = 100, and 1000 data points corresponds to m = 1000. Once you have formulated A and b, you can find the least squares solution by minimizing norm(A*x-b) in CVX, or you can find the least squares solution x using the baseline MATLAB as

A\b

or

pinv(A)*b

CVX will “come into its own” for this problem if you are adding, in a CVX (DCP) compliant manner, extra term(s) to the objective function and/or constraint(s) on x.

3

Thank you for your reply Mark! I was able to get this to work. In addition, I was able to implement a cost function to reduce the likelihood of overfitting by adding a second norm term.

Instead of using minimize( norm(Ax-b)), I used minimize( norm(Ax-b)+eta*norm(x) ) where \eta is an integer to help determine the cost function.