14

If what you want is to solve $\min \| Ax - b \|_{2}^{2} + \lambda^{2} \| x \|_{2}^{2}$ subject to $x \geq 0$, then this is easily implemented. Construct a matrix $C=\left[ \begin{array}{c} A \\ \lambda I \end{array} \right]$ and a vector $d=\left[ \begin{array}{c} b \\ 0 \end{array} \right]$. Then use your nonnegative least squares solver on $...


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

The standard tolerance for forming a pseudoinverse is to only invert singular values that are at least $\max(m,n) \epsilon \|A\|_2$, where $A$ is $m \times n$, $\epsilon$ is the machine precision, and $\|A\|_2$ coincides with the largest singular value of $A$. With that said, as J.M. mentioned, it is much more stable to avoid forming $A^H A$: First, we ...


11

Augh!! No, no, a thousand times, no! The reason people use SVD is precisely to avoid having to form the cross-product matrix $\mathbf A^\top\mathbf A$, since the formation of this matrix is a nice recipe for forming ill-conditioned linear systems! The decomposition is meant to be applied directly to $\mathbf A$. (See also some of my previous answers.) I ...


9

The Levenberg-Marquardt method can be used to minimize any problem of the form: $ \min f(x)=\sum_{i=1}^{m} f_{i}(x)^{2} $ However, if the objective functino to be minimized is not a sum of squares, then the method is no longer applicable.


8

$\mathbf{A}$ is an $(n+1) \times (n+1)$ matrix. It can be obtained as follows: $\textbf{A} = \left[ \matrix{1 & 1 & 1 & \cdots & 1 \cr x_0 & x_1 & x_2 & \cdots & x_{n} \cr \vdots & \vdots & \vdots & \cdots & \vdots \cr x_0^n & x_1^n & x_2^n & \cdots & ...


7

A few thoughts on your question: How you report your model fit will depend very much on your audience, and your field. For example, in my field, model fit statistics like R^2 are very rarely reported - regarded as neither impressive nor particularly useful. Instead, some criteria for how you arrived at the model you arrived at tends to be described, and ...


7

What you want to do here is partition N observations into K clusters who exhibit similar properties. This is called clustering and you can find more info here. Since you already have a numerical similarity measure, this makes me think about using the K-Means algorithm, in which you operate in several steps: Initialize cluster centroids randomly Assign each ...


6

The problem here is that the condition number for your monomial basis becomes very large, meaning that the value of the monomial is very sensitive to the value of the coefficient. Thus, when you try to compute in this basis for high degree polynomials, your answer is very inaccurate. There are polynomial bases that are better behaved, like any of the ...


6

First of all, rates of convergence are usually given in the form $$ \|u-u_h\| \leq C N^\alpha,$$ rather than equality. Furthermore, rates are asymptotic, i.e., only have to hold for $N\to \infty$. This means that you're unlikely to find a single $C$ and $\alpha$ such that your equation holds. Another reason why your approach doesn't work is that what you're ...


5

I colleague has found the Fast Hough Transform in the Gandalf library, which looks very promising but may be a lot of work to integrate, so I am looking for other approaches. The Gandalf implementation is interesting: they evaluation the accumulator space in a recursive way as if traversing a quad- or oct-tree. Regions without much density are thrown out as ...


5

I think what you are looking for is called "cluster analysis" or "clustering". Many different algorithms exist. In your case, you would want some "connectivity clustering", i.e. group elements together based on a property that links each two. Have a look at the clustering algorithms in scikits.learn (Python code) and the references mentioned there.


5

You're trying to fit a power law to your distribution. Very interesting. These show up all the time in graph theory, social networks, and a slew of other places. There's some tutorials on fitting your data here and here. Also, in reference to question A., how does the probability of a person buying land depend on how much land they already have? You may be ...


5

I wrote a Python package called PyPGE. PyPGE is a Symbolic Regression implementation based on Prioritized Grammar Enumeration (1), not Evolutionary or Genetic Programming. It produces a deterministic Symbolic Regression algorithm. (1) Worm, Tony, and Kenneth Chiu. "Prioritized grammar enumeration: symbolic regression by dynamic programming." Proceedings of ...


5

Numerical judgement of model choice: You model 36 observations with a model consisting of 12 or 13 predictor variables. This is most likely not a good model. Even if you reach a high $R^2_{adj}$, you most likely model a random pattern. Try to compare a computed $AIC$ (Akaike information criterion) or $BIC$ (Bayesian information criterion) of this model to ...


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

You could either fit a logistic function (possibly composing it with a linear function), use segmented regression, or classification and regression trees, among other options. The original data, shown in the figure below, was fitted in Gnuplot using the following commands: h(x) = k * 0.5 * (1.0 - tanh(0.5 * (a * x + b))) + c * x + d fit h(x) 'plot-EV.txt' ...


4

Based on your last comment (the fact that you might have multiple measurements and you don't care whether a certain set of measurements is above or below the fit), I think what you are looking for is a spline fit. You can do this using the scipy.interpolate B-spline routines. The script below generates three sets of data based on a function (that you can ...


4

You can easily do this with for loops. Just use n to define your loop limits. Here's a fully functional solution for MATLAB (just define n and x first): %pre-allocate A A = zeros(n+1); %first row: for j=0:n A(1,j+1)=sum(x.^j); end %rows 2 through n for i=1:n A(i+1,1:n)=A(i,2:n+1); %copy from previous row A(i+1,n+1)=sum(x.^(n+i)); %compute last ...


4

Let's use the example (approximate the square root function on $[0.25,1.0]$ with a quartic polynomial) to step through your calculations. I suspect that the code is going to work with only modest changes. Having chosen an initial set of six interpolation points: $$ x_k = 0.25, 0.4, 0.55, 0.7, 0.85, 1.0\;\; (k=0,\ldots,5) $$ we proceed to interpolate the ...


4

I was surprised for not receiving a satisfactory answer to the question above and my investigations showed me that this, indeed is an unexplored area. Hence, I put some effort in developing solutions to this problem, and published the following manuscripts: T. Birdal, B. Busam, N. Navab, S. Ilic and P. Sturm. "A Minimalist Approach to Type-Agnostic ...


3

If very different parameter sets give quite similar curves, you are most likely fitting too many parameters, which means that some of the fitted parameters will be sensitive to changes in the data. Look at the inverse Hessian of the objective function of the two fits to see how well the estimated parameters are determined by the data. Large diagonal ...


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

In case of an inequality constraint only and semidefinite $P$, your problem is convex, and there may be better alternatives (CVX mentioned in the answer by Victor Liu, or the methods of arXiv:1009.2065 (which has at thre end a reference to a an implementation). If $1<k<\infty$, you can use standard nonlinear programming software (see, e.g. http://neos-...


3

Without the equality constraint, the problem is convex, and any standard interior point convex optimization package can be applied to solve this efficiently, such as the high level modeling software CVX. With the equality constraint, the problem is no longer convex. However, you may approximate the solution by repeatedly solving a bunch of convex problems. ...


3

A general answer to the question is just impossible. This answer is limited to modeling the electrostatic term. We have recently published a detailed comparison of different methods to compute atomic charges, using two sets of 100+ penta-alanine conformers. (link) We tested both the robustness of the charges with respect to conformational changes and the ...


3

To ask about all ways is too ask too much, I guess. For the computationally most efficient way to encode the geometric information, see Section 2 of my survey paper Molecular modeling of proteins and mathematical prediction of protein structure, SIAM Rev. 39 (1997), 407-460.


3

By forcing your control points close to each other, your Vandermonde matrix becomes nearly singular and hence the condition number is poor and your solver fails. For simplicity, you have $x_2=x_1+\epsilon$, where $x_i$ are your control points. The two corresponding Vandermonde rows are: $$ \begin{align}\left(\begin{matrix} 1 & x_1 & x_1^2 &\...


3

Sections of Vandermonde matrices do tend to be ill conditioned. If indeed you were using the singular value decomposition for your fitting, you should have already seen a rash of tiny singular values. Did you remember to zero out those tiny singular values before computing the least squares solution? The ill-condition is easily visualized: if you plot ...


3

For the least-squares case anyway, one standard method is Linda Kaufman's variable projection algorithm (see this as well). This is based on earlier work by Gene Golub and Victor Pereyra for the separable nonlinear least squares problem, in which they show how the structure of the Moore-Penrose pseudoinverse can be exploited to yield algorithms that are ...


Only top voted, non community-wiki answers of a minimum length are eligible