18
The most important aspect of interpolation and curve fitting is to understand why high order polynomial fits can be an issue and what the other options are and then you can understand when they are/are not a good choice.
A few issues with high order polynomials:
Polynomials are naturally oscillatory functions. As the order of the polynomial increases, the ...
14
If a given nonlinear system is the first order optimality condition for an optimization problem, then we can often produce a more robust algorithm by using that information. For example, consider the equation
$$ f(x) = x^2 - \exp\big(-4(x-2)^2 \big) \qquad \text{[click for Wolfram Alpha]}$$
This clearly has a unique minimum and we expect our optimization ...
13
Inge Söderkvist (2009) has a nice write-up of solving the Rigid Body Movement Problem by singular value decomposition (SVD).
Suppose we are given 3D points $\{x_1,\ldots,x_n\}$ that after perturbation take positions $\{y_1,\ldots,y_n\}$ respectively. We seek a rigid "motion", i.e. a rotation $R$ and translation $d$ combined, applied to points $x_i$ that ...
11
I wrote a full answer (below the line) before discovering CVXPY, which (like CVX for MATLAB) does all the hard stuff for you and has a very short example almost identical to yours here. You only need to replace the relevant line with
p = program(minimize(norm2(A*x-b)),[equals(sum(x),1),geq(x,0)])
My old answer, doing it the harder way with CVXOPT:
...
10
Since one of my answers has been cited, I'll try to clarify why I suggested using IPOPT instead of MINPACK.
My objections to using MINPACK have nothing to do with the algorithms that MINPACK uses and everything to do with their implementation. My main objection is that the software dates back to 1980, and was last updated in 1999. Jorge Moré is retired; I ...
9
The problem is called Wahba's problem, one algorithm for it is called Kabsch algorithm, and the later more popular is called Davenport q method. It's apparently used and studied in aeronautics to determine a craft orientation. There are lots of reviews about the methods.
Beware that the best fit may include reflection.
Kabsch method computes a 3x3 ...
7
For solving the least squares problem in general as principal methods there are (matrix $A$ with full rank):
solve the system of normal equations $A^{T}Ax = A^{T}b$
use QR factorization
use SVD decomposition
Generally speaking the QR factorization is a method with a good balancing between accuracy and computational cost.
Using the normal equations is ...
6
In a word, yes.
A least squares problem arises in the context of an over determined system. That means that you have more equations than unknowns and would be the case if you had more data points than Fourier coefficients that you're trying to find. This would give you the "best fit" Fourier series of a given order. For example, you could find the best fit ...
6
David Eberly has a good description of the solution (with pseudocode) here: http://www.geometrictools.com/Documentation/IntersectionOfCylinders.pdf
The short summary: like most convex-convex collision detection algorithms, you search systematically for a separating axis between the cylinders.
For future reference on similar problems, the authors of Real-...
6
If you change variables to optimize for the residual of the linear part, then the Hessian will be a low-rank update to the identity. Then L-BFGS would work very well. Specifically, your problem takes the form
$$
\min_x \frac{1}{2}\|Ax-b\|^2 + \frac{\mu}{2}\|g(x)\|^2
$$
where $Ax=b$ is the linear PDE and $g$ is the nonlinear part, and $\mu$ is a tradeoff ...
5
I can't quite read your notation, but suppose we tried to find a set of parameters $x$ that give rise to predictions $y(x)=Ax$ by comparing with measurements $\bar y$. Then, if you try to find your parameters $x$ by solving
$$
\min_x \frac 12 \|\bar y - y(x)\|^2 = \frac 12 \|\bar y - Ax\|^2
$$
then it is indeed true that your covariance matrix is $(A^TA)^{-...
5
MATLAB has a built-in function lsqnonneg() which is an implementation of the active set method described in the book "Solving Least Squares Problems" by Lawson and Hanson (1974)
The \ solution of linear least squares problems in MATLAB is done by a variety of different algorithms depending on the exact structure of the problem. None of them are ...
5
In general, you can formulate this as a nonlinear least squares problem. If your values are known at points $(x_{i},y_{i})$, and the known values are $f_{i}$, then you can minimize
$\min_{a,b,c,d,e} \sum_{i=1}^{n} \left( f_{i} - a(bx_{i}+c)^{dy_{i}+e} \right)^{2}$
The Levenberg-Marquardt method is commonly used to solve nonlinear least squares problems ...
5
How would you solve the problem if you didn't need to do the centering?
Since $A$ is large and sparse, you'd probably pick an iterative method such as CGNE which depends on being able to perform matrix-vector multiplies $Ax$ and $A^{T}y$. It turns out that you can still use the same iterative method for the centered version of the problem since matrix-...
5
If $A$ is of full column rank and $A^{T}A$ is non-singular and well-conditioned, then you can compute the pseudoinverse as:
$ A^{\dagger}=(A^{T}A)^{-1}A^{T} $.
This will be faster than computing a QR or SVD factorization of $A$ but be careful about the conditioning of $A^{T}A$.
4
Instead of the problem you solved, solve
\begin{alignat}{1}
& \min_{x}\|Ax - b\|^{2}_{2}, \\
\mathrm{s.t.} & \quad\sum_{i}x_{i} = 1, \\
& \quad x_{i} \geq 0, \quad \forall{i}.
\end{alignat}
This problem is a differentiable, convex, nonlinear optimization problem that can be solved in CVXOPT, IPOPT, or any other convex optimization solver.
4
I'd guess a QR decomposition is better than solving the normal equations and faster than SVD.
There are some class notes that compare the three approaches.
Also: The QR decomposition for least-squares solution of overdetermined linear systems has been previously discussed here.
4
First, never calculate an explicit inverse, allthough sometimes you cannot avoid it. Most of the time however, the solution can be restated as multiplication of a vector with an inverse matrix. In that case several methods exist, which don't exlpicitely calculate the inverse.
That said, common regularization techniques are truncated singular value ...
4
If I read this correctly, you can rewrite your problem as an ordinary least square problem $\left\|Ax-b\right\|$, where A is the stacked Vandermonde-Matrix of the derivatives of the Zernike polynomials and $b$ the stacked measured data.
$$\chi^2 =\left \|\begin{pmatrix}
\frac{\partial P_n}{\partial x}\\
\frac{\partial P_n}{\partial y}
\end{pmatrix}C-\...
4
This looks isomorphic to Tikhonov regularization (a.k.a. ridge regression). Some googling brought up LSQR and the newer LSMR. Those links both have implementations in a number of languages. For large scale problems, Petsc has KSPLSQR built in.
Depending on what you mean by 'large scale', mlpack may also work. Mlpack has a tutorial covering both its ...
4
Pretty sure that what you want is Non-negative least squares. Plenty of implementations already exist, so you shouldn't have to write any code to get going. https://en.wikipedia.org/wiki/Non-negative_least_squares
4
@Stellos already gave the correct answer, but let me try to back it up with a bit of intuition:
Think of your function $f(x)$ as a function that would actually make sense for any real-valued argument $x$. Then, if that function happened to be of the form $f(x)=\sin(1000x)$, you would have lots of maxima and minima between each integer point and, in essence, ...
4
If we let
$\phi(x)=\sum_{i=1}^{m} F_{i}(x)^{2}$,
we could compute $\nabla \phi(x)$ by finite difference approximation. However, it is generally smart to make use of the special structure of $\phi(x)$ and derivatives of $F_{i}(x)$. Using the chain rule and exploiting the sum of squares structure of $\phi(x)$ we obtain that
$\nabla \phi(x)=2J(x)^{T}F(x)...
4
If you want to solve such problems with least squares, you need to first address two things:
You are interested in a continuous problem, not one particular discretization of it. Your statement "quite a fine mesh, the linear operator $\Delta$_h will be much larger than the 10 nonlinear equations at the end" shows that you are somehow aware of the problem, ...
4
$$\big\| \mathrm A \begin{bmatrix} \mathrm X & \mathrm X^2\end{bmatrix} - \begin{bmatrix} \mathrm B_1 & \mathrm B_2\end{bmatrix} \big\|_{\text{F}}^2 = \underbrace{\| \mathrm A \mathrm X - \mathrm B_1 \|_{\text{F}}^2}_{=: f_1 (\mathrm X)} + \underbrace{\| \mathrm A \mathrm X^2 - \mathrm B_2 \|_{\text{F}}^2}_{=: f_2 (\mathrm X)}$$
Everybody knows that
...
3
In the bound-constrained case, MATLAB defaults to using a trust-region reflective method found in Coleman, T.F. and Y. Li, "A Reflective Newton Method for Minimizing a Quadratic Function Subject to Bounds on Some of the Variables," SIAM Journal on Optimization, Vol. 6, Number 4, pp. 1040-1058, 1996. You might try that method, and see how it works out.
For ...
3
If you have access to the MATLAB optimization toolbox then this can easily be done using the quadprog() function. You'd start by writing the objective in quadratic form as
$ \| Ax - b \|_{2}^{2} = x^{T}(A^{T}A)x-2(A^{T}b)^{T}x+b^{T}b$
then multiply to get $P=A^{T}A$ and $q=-2A^{T}b$. Then your objective is
$f(x)=x^{T}Px+q^{T}x+b^{T}b$
and ready to ...
3
So, there are two ways to solve this problem.
The easy, non-rigorous way to solve this problem is to create a function that calls a MATLAB ODE solver using $K$ as a parameter, and returns the solution $(x_{1}(t), x_{2}(t))$ for all times corresponding to measured data points. Then use this function to construct a sum-of-squared error objective function, and ...
3
Sure, go with inverting $A^TA$. That matrix is symmetric and positive definite, so you can do that with CG as well. Alternatively, you can solve the original problem by computing a decomposition of the original, rectangular matrix, using an algorithm like QR.
3
If need be, you can use the "seminormal equations"; the key here is that the $\mathbf R$ factor in the QR decomposition of $\mathbf A$ is the same (up to changes in signs of entries) as the Cholesky triangle of $\mathbf A^\top\mathbf A$, since
$\mathbf A^\top\mathbf A=(\mathbf Q\mathbf R)^\top(\mathbf Q\mathbf R)=\mathbf R^\top\mathbf Q^\top\mathbf Q\mathbf ...
Only top voted, non community-wiki answers of a minimum length are eligible
