20

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 $...


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 & ...


8

My first guess was that they computed a minimax best-approximation polynomial to $\sqrt{x}$ on [0,5,1] with something like the Remez algorithm. However, plotting the difference $p(x)-\sqrt{x}$ I do not see 4 points of equioscillation, so this is not the polynomial that minimizes $$\max_{x\in [0.5,1]} |p(x)-\sqrt{x}|.$$ Maybe that 0.001 that they tried to ...


6

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 ...


6

What is the size of your $A$ matrix? Is $A$ sparse? Does $A$ have some other special structure? How many values of $\lambda$ do you want to try? Normally, you'd use the Cholesky factorization of $A^{T}A+\lambda I$ rather than the LU factorization since $A^{T}A+\lambda I$ is symmetric and positive definite. However, updating the Cholesky or LU ...


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 ...


6

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 ...


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

I once started writing anopen source version of Eureqa in Java. The project has limited capabilities but it implements the fitness function described in [1] and couple optimizations mentioned by the authors in other publications (e.g., searching for solutions in Pareto front). Link: https://github.com/pkoperek/hubert [1] Schmidt, Michael, and Hod Lipson. "...


4

After a cursory google search on the subject, it appears that "symbolic regression" is a problem that lends itself greatly to stochastic optimization algorithms like genetic programming (GP). It is conceivable that you should look for an open source genetic programming library with modules specifically for symbolic regression, such as DEAP (Distributed ...


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

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 ...


3

Given you are trying to find a path $f(t) = a + b t$ for each 3D component of an object to define its trajectory, you can formulate a Least Square problem to find the values for $a$ and $b$ based on $N$ pieces of data. The goal of the least square problem is to minimize the following cost function with respect to $a$ and $b$: $$J = \frac{1}{2} \sum_{i=1}^{...


3

Take a look at the book by Nocedal and Wright, "Numerical Optimization", to see the Levenberg-Marquardt method in more context than just fitting.


3

This is a constrained minimization problem (with linear constraints). You can write it like $$\min_X \frac{1}{2}\|AX-B\|^2_2 \quad \text{ s.t. } e^TX = 1$$ where $e$ is the vector of all ones. This problem is typically solved by forming the Lagrangian $L(X,\lambda)$ $$L(X,\lambda) = \frac{1}{2}\|AX-B\|^2_2 + \lambda(1-e^TX)$$ for which the constrained ...


3

David Ketcheson has already indicated the problem in his comment. I will flesh it out here. Note the form of the argument of the log in the logistic regression objective: $log[1+exp(−b_iA^T_ix)]$ where $b_i$ and $x$ are vectors and $A_i$ is a matrix. Only in the above order of matrix-vector multiplications will you get a scalar as an exponent. The Matlab ...


3

I found the gramEvol R package flexible and easy to use. They have a small tutorial in which they rederive Kepler's third law from data. Note that it relies on Genetic Programmic for its optimisation and thus might return different results if you run it twice.


3

An assortment of curves for fitting chemistry examples is presented in these Colby College class notes. Of particular application is the sigmoid response curve with variable "slope" for the central part of the curve: $$ f(x) = \frac{a}{1 + e^{bx - c} } + d $$ [This is similar to the suggested logistic function proposed in the first Answer, but has four ...


3

That might also have been a trick question. Let's say you want to solve the normal equations for $Ax=b$, i.e., $(A^T A) x = A^T b$. Let's assume for a moment that the questioner meant that $A$ is actually already stored in memory, so we know that that much memory is already available. Let's also assume that $A$ is tall and narrow (more specifically, has ...


3

If my computer history is not wrong, no one would use the initial guess obtained from the quadratic polynomial you provided and use the algorithm you suggested. This is due to two reasons: floating-point division was much more expensive than multiplication (up to 32x -or more if implemented in software- compared to 4x today) and obtaining the initial guess ...


3

This polynomial $p(x)$ solves the minimax optimization for $$\frac{p(x)}{\sqrt{x}}-1$$ over the interval $[\tfrac12,1]$.


2

Minimizing the 2-norm of $x$ among all least squares solutions is relatively easy to do- this is the pseudoinverse solution. It can be computed using either a rank revealing version of the QR factorization of $A$ (there are specialized versions of this that work well on large and sparse $A$ matrices) or by using the Singular Value Decomposition. Since you'...


2

There is also a package for R called rgp. Visit this link. https://cran.r-project.org/web/packages/rgp/index.html I have not used rgp as I have only begun to use R seriously but it seemed like a good lead. I have another one for you that looks really promising but I have a mac and cannot use it: http://dev.heuristiclab.com/wiki/AdditionalMaterial/ECML-...


2

What are good parametrizations of rational functions for response surface models? A widely used flexible parametrization of (piecewise) rational functions are non-uniform rational basis spline (NURBS) models. The advantage of this parametrization is that the effect of the parameters can be understood intuitively. The effect of the denominator here is to ...


2

If you know that you have $n$ points $X_1,...,X_n$ and the desired pairwise distances $d_{ij}$ between $X_i$ and $X_j$ you could try and optimize the functional $$ \sum_{i<j} (d_{ij}-|X_i-X_j|)^2$$ In my opinion this looks good as an optimization problem. You can compute the gradient explicitly (even the Hessian, if you are motivated) and then use some ...


2

So the correct way is to use the probability distribution of your $c$ in the modeling process. Imho, the most natural way of describing and solving this type of problem is the Bayesian approach $$P(f|X,m)=\frac{P(X|f,m)P(f|m)}{P(X|m)}.$$ Noise free $c$ You have the data points $X_i=(x_i,d_i)$. The model $m$ predicts the values of $d$: $$d(x)=m(x,f)+\...


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