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
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 ...
10
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
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
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
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
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)...
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
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 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
...
4
As the commenters have already said, for your specific problem, QR factorization might be a better approach not because of speed but because of numerical stability.
It's still a good question to ask in general.
One of the advantages you cite is that $LDL^*$ can be used for indefinite matrices, which is definitely a point in its favor.
The linear algebra ...
4
SYRK is not really relevant here I think; it is just something else that happens to have the same name "rank k update".
In your case, you need to know how to update a QR factorization by inserting rows; a good reference is Golub, Van Loan, section 6.5.3: Appending or Deleting a Row.
Many computational environments have it already implemented for you, see e....
4
If you have a good LP solver, then the linear programming approach often works well. You don’t want to implement your own simplex or interior point code for LP though.
A specialized variant of the simplex method due to Barrodale and Roberts is a popular approach to this problem, but it takes some effort to implement this efficiently and you might be better ...
3
Hint:
If the two cylinders are parallel, the problem is easy. Otherwise, if the perpendicular distance between the axis exceeds the sum of the radii, there is no intersection.
Otherwise, there is an intersection curve between the two lateral surfaces. It can be expressed analytically: without loss of generality, the second cylinder is vertical with a base ...
3
The set of convex sets is infinite dimensional, which is not amenable to computation. I would discretize this by replacing your solution set (consisting of all convex sets) by a finite dimensional set -- for example, polygons with $N$ vertices. This finite dimensional set can conveniently be represented as the intersection of $N$ half-spaces and you would ...
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
I'll add a real response here, following the comment I just left. As @GeoffOxberry suggests, you might be able to use active subspaces to, in essence, preprocess your objective function and eliminate (linear combinations of) variables. Try the following first. Randomly sample your variables according to some density. One reasonable choice would be to sample ...
3
You could try my MATLAB solver PDCO which uses an interior method and will be happy that your n < m.
Use
options.Method = 1 % the default
d1 = 1e-4
d2 = 1
3
You should take a look at the excellent book by Nocedal and Wright on "Numerical Optimization". It has a great deal of material on these kinds of problems.
The algorithm you refer to, where one only adds or subtracts one variable to the active set at a time is described in detail in the book. The reason for only adding or subtracting one variable at a time ...
Only top voted, non community-wiki answers of a minimum length are eligible
