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:
-
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.
-
Would I use the sum_square(A*x-b) function in cvx?
Thank you in advance for your help!
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.
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.
