# Least Square Method with Positive Parameters

this is my first post here in the Stack Exchange. A friend told me about this forum and I'm giving it a try.

I searched a bit past threads, but couldn't find what I wanted, so I'm posting the problem here.

So, what I have is the following:

I have an $n(x) = \frac{e^{-x^2/2}}{\sqrt{2\pi}}$ in $\mathbb{R}$ and for an Smoothed Particle Hydrodynamics simulation I need to find a set of $x_i$ and $\nu_i$ such that I have

$$n_{SPH} = \sum_{i=0}^N \nu_i W(x-x_i,h)$$

The closest I can to the function $n(x)$. The $W(x-x_i,h)$ is a dirac kernel, i.e., (in principle) Compact Support, $\int_\mathbb{R} W(x,h) dx= 1$ and $\int_\mathbb R W(x-x',h)f(x')dx' \rightarrow f(x)$ as $h\rightarrow 0$. This could be an gaussian, a spline or something like that. Also, this kernel is positive and simetric.

Idealy, I would do Least Squares in both $x_i$ and $\nu_i$, but lets just say that I alredy have the $x_i$. I would have to miniminize the function

$$f(\nu_1,\nu_2,\ldots,\nu_N) = \int_\mathbb R dx (n(x)-n_{SPH}(x))^2 =\int_\mathbb R dx (n(x)-\sum_{i=0}^N \nu_i W(x-x_i,h))^2$$

Well, I wrote my program to do the least squares according to this. The problem is that I have an aditional constraint that is to have $\nu_i \geq 0$, and the coeficients that it spits aren't all positive, and this is causing me problems.

I don't know much about convex algorithms and such that I found in other threads. I'm an undergraduate physics my knowledge is restricted basicaly to C(not C++, not python, not fortran) progaming, a bit of numerical calculus and some other simple algorithms. About libraries, I know how to use C's Standard Library and (can handle) Gnu Scientific Library, but I'm not realy familiarized with lots of numerical packages.

If possible, I wanted a reference or an explanation how to code a solution to this, I have no preference nor dislike to any ready solution.

Edit - 28/12/2011:

Thanks a lot J.M., sorry for the delay in answering, but I have bumped in another problem.

In order to use the algorithm listed on the link, I have to cast it on the form required by the book.

problem stated in the link is to miniminize $||Ex-f||$ with the constraint that $x\geq 0$.

In my case I have something like miniminizing $g(\vec x) = N_0 - <\vec N,\vec x> + <\vec x, M \vec x>$

where

$$M_{ij}=\int_\mathbb R W(x-x_i,h)W(x-x_j,h)$$ $$N_i = \int_\mathbb R n(x)W(x-x_i,h)$$ $$N_0 = \int_\mathbb R n(x)^2$$

So, casting in the problem form, there would be the following identifications: $$M = E^T E$$ $$\vec N = 2 E^T f$$ $$N_0 = ||f||^2$$

So my question is the following:

Is there any way to assure that the matrix E is uniquely determined, and if so, a simple way to calculate it.

I was thinking on cholesky decomposition, but I'm not sure if it would work.

-
You'll want to see this. – J. M. Dec 11 '11 at 2:47
See also this and this. – J. M. Dec 11 '11 at 2:51
J. M. I did what I just said, I tried to use cholesky decomposition on the M martix, and I couldn't do it because it's not positive definite, GSL returned error because of that. Any sugestion as how to find that E matrix? – Hydro Guy Jan 5 '12 at 20:13
@Spirit_of_the_Air Do you know any Matlab? If so you'll want to check this solver: stanford.edu/group/SOL/software/pdco.html. Regarding the Cholesky decomposition of $M = E^T E$, it may be obtained from the $QR$ decomposition of $E = QR$ where $Q$ has orthogonal columns and $R$ is upper triangular, since then $M = E^T E = R^T Q^T Q R = R^T R$. Since the Cholesky decomposition is unique, this is it. – Dominique May 25 '12 at 21:42