# Constrained Least Squares

1. Aug 23, 2007

### hotvette

Anybody know the math/theory behind linear least squares where the curve is forced to go through the first and last data points?I'm specifically dealing with cubic polynomials.

In standard linear least squares formulation (i.e. ATAc = ATy) the curve doesn't, in general, go through any of the data points. I've actually managed to come up with a formulation based on Bezier math that allows standard least squares treatment to be used, but the development is problem dependent and rather long and tedious.

Just wondering if there is an easiser way. Ideas?

Last edited: Aug 24, 2007
2. Aug 23, 2007

### CRGreathouse

A cubic polynomial is of the form y = (x - a)(x - b)(x - c)d. If you want it to go through two particular points, that constrains two of the variables -- just adjust the others as needed.

3. Aug 24, 2007

### hotvette

Hmmm, not sure I follow. Picking a and b will set 2 of the zeros of the cubic, meaning it goes though the points (a,0) and (b,0), but I'm trying to force it to go through (x0, y0) and (x1, y1).

After a-l-o-t of work I was able to come up with:

y = (c(3x3-27x2+72x-48) + d(-3x3+18x2-27x+12) + x3+6x2-42x+62)/27

for (x0, y0) = (1,1) and (x1, y1) = (4,2).

Let me re-word the question. Is there a relatively straightforward way to construct an equation of a cubic polynomial that goes though 2 specific points and has 2 parameters to adjust? I need it in this form to evaluate ∂y/∂c and ∂y/∂d for least squares treatment. Unless there is another way..........

Last edited: Aug 24, 2007
4. Aug 24, 2007

### AlephZero

Yes.... but if two of the constants a, b, c are complex, then it's not a very good place to start from. For example, x^3 + x gives a = 0, b = i, c = -i.

Bezier or Bernstein polymonials (apologies, can't remember which is the right name without looking it up) seem a better idea.

If your cubic goes through $$(x_0,y_0)$$ and $$(x_1,y_1)$$

$$y = a(x-x_0)^3 + b(x-x_0)^2(x-x_1) + c(x-x_0)(x-x_1)^2 + d(x-x_1)^3$$

$$y_0$$ and $$y_1$$ fix the values of a and d. Find b and c by least squares.

5. Aug 24, 2007

### hotvette

Perfect, thanks! I somehow frequently manage to come up with the most complicated way to do things...

$$y = a(x-1)^2(x-4) + b(x-1)(x-4)^2 + \frac{(x+2)}{3}$$

is much easier.
-----------------

Question - where'd you find that form of a cubic poly?

Last edited: Aug 24, 2007
6. Aug 24, 2007

### AlephZero

Try Googling for the B-splines or Bernstein polynomials, (the B in B-splines probably stands for "Bernstein"). Bezier splines are very closely related. The "slope" control point at one end of a Bezier spline depends on a and b, and at the other on c and d.

I first saw this in some tutorial notes on splines, long before the Internet existed - I can't remember exactly where.

7. Aug 27, 2007

### nezard

I am sitting on the same problem for several day already and follow this discussion with the great interest :)

Now, if you fix values of a and d (with previous statement), how do you apply/formulate LSS for remaining points?

Thx

8. Aug 27, 2007

### hotvette

Welcome to the forums.

There may be multiple ways to approach it, but here is what I did.

1. Pick arbitrary values for b and c
2. Based on the values of a,b,c,d calculate f(x) and (f-y) for each point
3. Form J by calculating ∂y/∂b and ∂y/∂c for each point
4. Solve JTJs = JT(f-y) for s
5. Update values for b & c: [b' c'] = [b c] - s

It only takes 1 iteration because the problem is linear. Net is to apply non-linear LS technique. Check out the attachment in the last post of this thread:

.

Last edited: Aug 27, 2007
9. Aug 28, 2007

### nezard

To be honest, I also do not see how did you evolve here

$$y = a(x-1)^2(x-4) + b(x-1)(x-4)^2 + \frac{(x+2)}{3}$$

$$y = a(x-x_0)^3 + b(x-x_0)^2(x-x_1) + c(x-x_0)(x-x_1)^2 + d(x-x_1)^3$$

BTW, it would be more then wonderful if you could give us parts of your example that you did on Bezier CPs.

many thx

10. Aug 29, 2007

### hotvette

Yeah, it's a bit convoluted. What I wanted was a form where y = y1 at x = x1 and y = y0 at x = x0. If you evaluate the expression at x = x0 and x = x1 you get:

$$y_1 = a(x_1-x_0)^3 \ \ and\ \ y_0 = d(x_0-x_1)^3$$

which means

$$a = \frac{y_1}{(x_1-x_0)^3} \ \ and\ \ d = \frac{y_0}{(x_0-x_1)^3}$$

Substituting those expressions:

$$y = y_1\frac{(x-x_0)^3}{(x_1-x_0)^3} + b(x-x_0)^2(x-x_1) + c(x-x_0)(x-x_1)^2 + y_0\frac{(x-x_1)^3}{(x_0-x_1)^3}$$

I then discovered (thanks to HallsofIvy) that the cubes aren't really needed on the 1st and last terms:

$$y = y_1\frac{(x-x_0)}{(x_1-x_0)} + b(x-x_0)^2(x-x_1) + c(x-x_0)(x-x_1)^2 + y_0\frac{(x-x_1)}{(x_0-x_1)}$$

Substituting:

$$x_0 = 1\ \ \ y_0 = 1\ \ \ x_1 = 4\ \ \ y_1 = 2$$

and simplifying will yield the final form.

Re your second question, what are Bezier CPs? Cubic Polynomials? What specifically would you like to know?

Last edited: Aug 29, 2007
11. Aug 31, 2007

### nezard

However, in the meantime I managed to discover another way to approximate “my” curve, yielding less computation overhead if compared to the least squares. To put it straight, my input is not discrete set of data (points), but the full mathematical expression y=y(x), or even parametric y=y(t) & x=x(t). The goal of an approximation here was both performance and visual-quality driven, e.g. to avoid splitting the curve into arbitrary number of very short segments interconnected with the straight line, giving overall appearance of one smooth curve, where each point is calculated by a dedicated call to above stated functions. Instead, I wanted limited number of fine cubic splines, or more precisely Bezier cubic splines, to approximate everything. Bezier simply therefore, it is an integral part of java native graphics library (there is no other curve built in java) and it comes with the neat built in antialiasing mechanism for smoothness.

The problem unexpectedly arose by defining control points for Bezier, since they are not on the curve itself. As I was trying to find solution I stumbled upon this forum and several of your posts related to this particular problem :)

BTW, the CP actually meant Control Point for a cubic Bezier :)

12. Sep 7, 2007

### hotvette

OK, sounds like your task is to find Bezier control points to represent certain functions, like y = ex or y = sin(x), as accurately as possible.

Have you figured out how to do it, or do you need help? If you need help, tell is some specific function examples you are working with.

13. Sep 8, 2007

### nezard

Well, I do have a solution here, but as someone said “there is always room for improvement”

I'll describe briefly both, problem and the actual solution, and you are welcome to comment on it. Although I am not going into details, they can always be provided if needed. Also in case you would choose quite different approach, please feel invited to present it

(1) Assume function R = a + b cos(t) in polar coordinates, or in parametric form
x = a cos(t) + b cos(t)^2, y = a sin(t) + b cos(t)sin(t), which you want to approximate as close as possible with Bezier, subdivided in N sections, where N is usually in range 8-16 (to get one smooth and continuous curve).

(2) Obviously, every point on the curve is known, so start-point and end-point for Bezier is not an issue, regardless of how the overall range t=0-2Pi is subdivided into sections. The only problem here is how to find control points

(3) As the curve function is known, its derivative is known as well, which suits well the basic idea of Hermite approximation. A typical Hermite is defined with four points, having control points defined as the first derivative in its start, and in its end point, or more precisely B = {P1, P2, P1', P2'} for any arbitrary segment.

(4) Once having the Hermite approximation Points it is not a problem to transform them into Bezier of four points.

(5) The problem with Hermite approach however, is approximation accuracy, for the curve (segment) has good shape, but it still does not fit the original curve exactly, except in start and end points. Here I do the trick, by calculating the needed scale factor, required to multiply derivatives, for the curve to pass through one arbitrary point on the segment. As a result, we are getting acceptable approximation through three points per segment. There is of course still some minor remaining deviation between the original shape and the approximated one, but its the one I can well live with.

14. Sep 9, 2007

### hotvette

Your method sounds perfectly valid. Whether or not there's a 'better' approach depends on a number of things:

- whether you are enforcing C2 continuity at the connection points
- whether the calcs are done on the fly and therefore need to be computationally efficient
- I'm sure there are more

Just for kicks, I used a=1 and b=3. Also, because of symmetry, I only looked at t = 0 to pi. I broke the function into 4 segments, which is plenty. I used points where the function has vertical or horzontal slope as the Bezier connect points. After a fair amount of playing around, I came up with the following control points.

segment 1:
P0 = (4, 0)
P1 = (4, 1.234536617)
P2 = (3.013679806, 2.24436971)
P3 = (1.814333489, 2.24436971)

segment 2:
P0 = (1.814333489, 2.24436971)
P1 = (0.614987172, 2.24436971)
P2 = (-0.083333333, 1.314059632)
P3 = (-0.083333333, 0.493006649)

segment 3:
P0 =(-0.083333333, 0.493006649)
P1=(-0.083333333, -0.328046334)
P2 = (0.534044209, -0.840154844)
P3 = (1.102333177, -0.840154844)

segment 4:
P0 = (1.102333177, -0.840154844)
P1 = (1.670622146, -0.840154844)
P2 = (2, -0.406642476)
P3 = (2, 0)

Visually, I cannot tell the different between the connected Bezier segments and the function. Plot it and see for yourself. I'll tell you the details of what I did tomorrow. It's late and I'm tired.
:zzz:

Last edited: Sep 9, 2007
15. Sep 9, 2007

### nezard

For illustration only, I enclosed a comparison of your approximation and the real function of Limacon of Pascal, showing a fairly good match. The approximation, in magenta, also shows all Bezier points. This (a=1 and b=3) however, is not very difficult shape. The real “fun” begins when a=4, 5, 6 or 7 having b=3. Then you might “wish” to have more then four segments in range 0-Pi.

As far as conditions are considered, yes, it should be calculated on the fly (many times in a loop per second) and therefore the performance is an issue. And for the other, I haven't tried forcing the continuity at connection points, for I consider each segment separately, but somehow I imagine having difficulties with it on the span of four segments anyway.

I am looking forward seeing the details of your method

#### Attached Files:

• ###### Approx1.zip
File size:
124.2 KB
Views:
57
16. Sep 9, 2007

### hotvette

Since computation performance is important, I don't think there is anything I can suggest that is better than what you are currently doing. What I did was computationally intensive and doesn't appear to have resulted in anything better. I'm surprised at the accuracy you are trying to achieve if the example I gave is considered a 'fairly good match'.

Anyway, here is what I did for a=1 and b=3:

1. Split the curve for t=0 to pi into 4 segments.

2. The boundary of each segment is a point at which the slope is either horizontal or vertical.

3. Enforce C1 and C2 continuity at segment intersections. This results in 2 control points to adjust in the first segment and 1 control point to adjust in each of the remaining segments

4. Sum perpendicular distances at 10 equally distributed points in each segment

5. Using a simplex optimization program, adjust the control points in each successive segment by minimizing the qty calculated in 4.

It sounds like what you are doing is enforcing C1 (and C0 obviously) continuity at each segment boundary, which results in 2 control points to adjust for each segment. The control points are chosen to force an exact match at some point within the segment (probably around t = 0.5). I've seen articles on fitting a Bezier segment to a quarter circle that use the exact same approach. One article insists that the method is the only correct one.........

One thing I'm curious about is what method you are using to judge quality of fit. I can think of a least 4 ways. In the case of the 1st segment (a=1, b=3), I came up with 4 different sets of control points depending on the quality of fit criteria:

Px0 = Px1 = 4
Px3 = 1.814333489

Py0 = 0
Py2 = Py3 = 2.24436971

Case 1:
Px2 = 3.13839852
Py1 = 1.110593575

Case 2:
Px2 = 3.00408392
Py1 = 1.244587413

Case 3:
Px2 = 3.013679808
Py1 = 1.234536616

Case 4:
Px2 = 2.931618476
Py1 = 1.319853858

I'm curious as to which you consider a better fit and why.

Last edited: Sep 9, 2007
17. Sep 9, 2007

### nezard

Well, I did not mean to sound critical with 'fairly good match', but I'll rephrase my statement gladly into: For most practical use I can imagine it is indeed an excellent match. However, only if a < b, like a=1 and b=3. Since this is not always the case I enclose a snapshot of the curve when parameters are a=5 and b=3, which illustrate rather difficult shape to be approximated in the second segment, as well as dissatisfying approximation - with my method.

To be honest, seeing how they approximate 1/4 of the circle in java brought me to the idea to try doing the same with Limacon. Unfortunately, seeing it coded in java is not enough to be able to grasp an underlaying logic which makes it work the right way – in particular if you do not have strong or fresh math background like in my case. Luckily, I came across an excellent book of prof. David Salomon from CSU, which helped me gaining essential insight into the area.

Although not aware of articles you mention, yes I use point t=0.5 to adjust control points....

It was a pleasure of 'writing' to you - thx

#### Attached Files:

• ###### Approx2.zip
File size:
74 KB
Views:
53
18. Sep 9, 2007

### nezard

OK, from whatever reason I received only one half of you post at first and missed the rest, e.g. the question.

For this purpose there is no particular 'quality assurance' methodology I can think of. Visual observation and subjective judgment is the best I can offer at this point. Even though I would probably go for 3 or 4, all four cases would do for constant a=1 and b=3, as the misalignment does not seem to be that considerable (more then 1-2 pixel on the screen length) and all curves overall look smooth. Insofar I am curious about how do you evaluate quality of fit. BTW, have a look at the curve I am getting for four segments:

Segment 1:
P0 {4.,0.}
P1 {4.,1.26725}
P2 {2.98297,2.24437}
P3 {1.81433,2.24437}
Segment 2:
P0 {1.81433,2.24437}
P1 {0.695769,2.24437}
P2 {-0.0833333,1.38385}
P3 {-0.0833333,0.493007}
Segment 3:
P0 {-0.0833333,0.493007}
P1 {-0.0833333,-0.302675}
P2 {0.504082,-0.840155}
P3 {1.10233,-0.840155}
Segment 4:
P0 {1.10233,-0.840155}
P1 {1.62068,-0.840155}
P2 {2.,-0.451426}
P3 {2.,0.}

But the problem as I see it is not to choose from these four. It is rather to ensure the method to generate control points does at least as good results for other values of a and b. In other words, that you can relay on the method always to provide you with the 'fairly good match', in all ranges of a/b. If you have an algorithm for 1 and 3, give it a try with a=5 for instance.

19. Sep 11, 2007

### hotvette

OK, I think I understand what you are trying to do. Seems like quite a challenge to develop an algorithm that will be adequate for all possible combinations of a and b. I need to ponder this a bit.

Question - what technique are you using to adjust the control point values to force the curve to go though a mid point? Some sort of iterative method?

20. Sep 12, 2007

### nezard

The technique to adjust CPs is not iterative. An iteration would be to expensive I expect.

As I do not know how to edit formulas here, I enclosed a brief description of the method I use in pdf file.

File size:
43.7 KB
Views:
102