On Tue, Mar 5, 2013 at 6:23 AM, Charles R Harris <[email protected]>wrote:
> > > On Tue, Mar 5, 2013 at 12:41 AM, Jaime Fernández del Río < > [email protected]> wrote: > >> On Mon, Mar 4, 2013 at 8:37 PM, Charles R Harris < >> [email protected]> wrote: >> >>> >>> There are actually seven versions of polynomial fit, two for the usual >>> polynomial basis, and one each for Legendre, Chebyshev, Hermite, Hermite_e, >>> and Laguerre ;) >>> >> >> Correct me if I am wrong, but the fitted function is the same regardless >> of the polynomial basis used. I don't know if there can be numerical >> stability issues, but chebfit(x, y, n) returns the same as >> poly2cheb(polyfit(x, y, n)). >> >> In any case, with all the already existing support for these special >> polynomials, it wouldn't be too hard to set the problem up to calculate the >> right coefficients directly for each case. >> >> >>> How do you propose to implement it? I think Lagrange multipliers is >>> overkill, I'd rather see using the weights (approximate) or change of >>> variable -- a permutation in this case -- followed by qr and lstsq. >>> >> >> The weights method is already in place, but I find it rather inelegant >> and unsatisfactory as a solution to this problem. But if it is deemed >> sufficient, then there is of course no need to go any further. >> >> I hadn't thought of any other way than using Lagrange multipliers, but >> looking at it in more detail, I am not sure it will be possible to >> formulate it in a manner that can be fed to lstsq, as polyfit does today. >> And if it can't, it probably wouldn't make much sense to have two different >> methods which cannot produce the same full output running under the same >> hood. >> >> I can't figure out your "change of variable" method from the succinct >> description, could you elaborate a little more? >> > > I think the place to add this is to lstsq as linear constraints. That is, > the coefficients must satisfy B * c = y_c for some set of equations B. In > the polynomial case the rows of B would be the powers of x at the points > you want to constrain. Then do an svd on B, B = u * d * v. Apply v to the > design matrix of the unconstrained points A' = A * v.T so that B' becomes > u * d. The coefficients are now replaced by new variables c' with the > contraints in the first two columns. If there are, say, 2 constraints, u * > d will be 2x2. Solve that equation for the first two constraints then > multiply the first two columns of the design matrix A' by the result and > put them on the rhs, i.e., > > y = y - A'[:, :2] * c'[:2] > > then solve the usual l least squares thing with > > A[:, 2:] * c'[2:] = y > > to get the rest of the transformed coefficients c'. Put the coefficients > altogether and multiply with v^T to get > > c = v^T * c' > > There are a few missing `'` in there, but I think you can get the idea, we are making the substitution c = v^T * c'. Chuck
_______________________________________________ NumPy-Discussion mailing list [email protected] http://mail.scipy.org/mailman/listinfo/numpy-discussion
