Least Squares Problems with Absolute Quadratic Constraints

This paper analyzes linear least squares problems with absolute quadratic constraints. We develop a generalized theory following Bookstein’s conic-ﬁtting and Fitzgibbon’s direct ellipse-speciﬁc ﬁtting. Under simple preconditions, it can be shown that a minimum always exists and can be determined by a generalized eigenvalue problem. This problem is numerically reduced to an eigenvalue problem by multiplications of Givens’ rotations. Finally, four applications of this approach are presented.


Introduction
The least squares methods cover a wide range of applications in signal processing and system identification 1-5 . Many technical applications need robust and fast algorithms for fitting ellipses to given points in the plane. In the past, effective methods were Bookstein's conic-fitting or Fitzgibbon's direct ellipse-specific fitting, where an algebaic distance with a quadratic constraint is minimized 6, 7 . In this paper, we develop an extended theory of minimization of least squares with a quadratic constraint based on the ideas of Bookstein and Fitzgibbon. Thereby, we show the existence of a minimal solution and present the uniqueness regarding to the smallest positive generalized eigenvalue. So, arbitrary conic fitting problems with quadratic constraints are possible.
Let A ∈ R n×m be matrix with n ≥ m ≥ 2, C ∈ R m×m be symmetric matrix, and d ∈ R a real value. We consider the problem of finding a vector x ∈ R m which minimizes the function F : R m → R defined by The side condition x t Cx d introduces an absolute quadratic constraint. The problem 1.1 is not a special case of Gander's optimization as presented in 8 , because in our case C is a real symmetric matrix in contrast to the approach of Gander, where the side condition considers real symmetric matrices C t C which are positive definite. For our considerations, we require the following three assumptions Assumption 1.1. By replacing C by −C , we consider d ≥ 0. For d 0, the trivial solution x 0 ∈ R m fulfills 1.1 . Therefore, we demand d > 0.
Assumption 1.2. The set N : {x ∈ R m , x t Cx d} is not empty, that is, the matrix C has not been less than one positive eigenvalue. If we assume that C has only nonpositive eigenvalues, it would be negative semidefinite and holds. With d > 0 it follows that the set N would be empty in this case. Assumption 1.3. In the following, we set S A t A and assume that S is regular. S is sometimes called scatter matrix.
In the following two sections, we introduce the theoretical basics of this optimization. The main result is the solution of a generalized eigenvalue problem. Afterwards, we reduce this system numerically to an eigenvalue problem. In Section 5, we present four typical applications for conic fitting problems with quadratic constraints. These approximations contain the ellipse fitting of Fitzgibbon, the hyperbola fitting of O'Leary, the conic fitting of Bookstein, and an optical application of shrinked aspheres 6, 7, 9, 10 .

Existence of a Minimal Solution
Theorem 2.1. If S ∈ R m×m is regular, then there exists a global minimum to the problem 1.1 .
Proof. The real regular matrix S A t A is symmetric and positive definite. Therefore, a Cholesky decomposition S R t R exists with a regular upper triangular matrix R ∈ R m×m . In 1.1 , we are looking for a solution x ∈ R m minimizing F x Ax 2 x t Sx x t R t Rx subject to x t Cx d.

2.1
With R regular we substitute x by R −1 y for y ∈ R m . Thus, we obtain an equivalent problem to 2.1 , where we want to find a vector y ∈ R m , minimizing F y y t y y 2 subject to y t R −1 t CR −1 y d.

2.2
Now, we define G : R m → R with G y y t R −1 t CR −1 y − d and look for a solution y on the zero-set N G of G with minimal distance to the point of origin. Let y 0 ∈ N G / ∅ and K m r 0 0 be Journal of Applied Mathematics 3 the closed sphere of R m in 0 with radius r 0 y 0 . Because of y 0 ∈ N G ∩ K m r 0 0 and G being continuous, the set is nonempty, closed and bounded. Therefore, for the continuous function F y y 2 exists a minimal value y M on N G ∩ K m r 0 0 with For all y from N G \ K m r 0 0 , it is F y > r 2 0 . So, y M is a minimal value of F in N G . By the equivalence of 2.1 , and 2.2 the assumption follows.

Generalized Eigenvalue Problem
The minimization problem in 1.1 induces a generalized eigenvalue problem. The following theorem is already proven by Bookstein and Fitzgibbon for the special case of ellipse-fitting 6, 7 . Proof. Let G : R m → R be a defined as G x : d − x t Cx. For G x 0 and d > 0 follows x / 0. Further, G is continuously differentiable with dG/dx −2Cx / 0 for all x of the zero-set of G. So, if x s is a local extremum of F x subject to G x 0, then it is rank dG/dx x s 1. Since F is also a continuously differentiable function in R m with m > 1, it follows by using a Lagrange multiplier 11 : if x s is a local extremum of F x subject to G x 0, then a λ 0 ∈ R exists, such that the Lagrange function φ : R m 1 → R given as has a critical point in x s , λ 0 . Therefore, x s fulfills necessarily the equations: Journal of Applied Mathematics The first equation describes a generalized eigenvalue problem with With d > 0, x s / 0 and x s fulfills 3.6 , λ must be a generalized eigenvalue, and x s is a corresponding eigenvector to λ of 3.6 , so that S − λC is a singular matrix. If λ 0 is an eigenvalue and x 0 / 0 a corresponding eigenvector to λ 0 of 3.6 , then every vector αx 0 is also a solution of 3.6 for λ 0 . Now, we are looking for α, such that x s αx 0 satisfies 3.5 . For λ 0 / 0 and 3.4 , it follows Because the left side and the numerator are positive, the denominator must also be chosen positive, that is, only positive eigenvalues solve 3.4 and 3.5 . By the multiplication with λ 0 , follows and x s α · x 0 fulfills the constraint G x s 0.
Remark 3.2. Let x 0 be a generalized eigenvector to a positive eigenvalue λ 0 of problem 3.6 . Then are solutions of 3.8 .

Lemma 3.3. If S is regular and C is symmetric, then all eigenvalues of 3.1 are real-valued and different to zero.
Proof. With det S / 0, λ 0 / 0 in 3.1 . The Cholesky decomposition S R t R with a regular upper triangular matrix R yields 3.1 With R invertible and the substitution μ 0 λ −1 0 , y s Rx s , we obtain an eigenvalue problem to the matrix R t −1 CR −1 : Furthermore, we have Journal of Applied Mathematics 5 Therefore, the matrix R t −1 CR −1 is symmetric and all eigenvalues μ 0 are real. With λ 0 μ −1 0 for μ 0 / 0 follows the proposition.
Remark 3.4. Because of S regular and λ / 0 in we can consider the equivalent problem with μ 0 λ −1 0 instead of 3.11 : This system is called inverse eigenvalue problem to 3.1 . Here, the eigenspaces to the generalized eigenvalue λ 0 in 3.1 and to the eigenvalue 1/λ 0 in 3.14 are identical. Therefore, the generalized eigenvectors in 3.1 are perpendicular.
Definition 3.5. The set of all eigenvalues of a matrix C is called spectrum σ C . We call the set of all eigenvalues to the generalized eigenvalue problem in 3.1 also spectrum and denote σ S, C . σ S, C is defined as the set of all positive values in σ S, C .
Remark 3.6. In case of rg C < m rg S , the inverse problem in 3.14 has eigenvalue 0 with multiplicity rg S − rg C . Otherwise, for μ / 0 and μ ∈ σ S −1 C follows 1/μ ∈ σ S, C .
The following lemma is a modified result of Fitzgibbon 7 . Proof. With S being nonsingular, every generalized eigenvalue λ 0 of 3.1 is not zero. Therefore, it follows for the equivalent problem 3.11 that μ 0 λ −1 0 is also a positive eigenvalue to R −1 t CR −1 , where R is an upper triangular matrix to the Cholesky decomposition of S. With Sylvester's Inertia Law, we know that the signs of eigenvalues of the matrices R −1 t CR −1 are the same as those of C.
For the following proofs, we need the lemma of Lagrange see, e.g., 12 .
Then x s is a minimal solution of f in N g .
Definition 3.9. Let λ * be the smallest positive value of σ S, C and x * s a corresponding generalized eigenvector to λ * to the constraint x * t s Cx * s d.

6
Journal of Applied Mathematics Lemma 3.10. Let S − λC be a positive semidefinite matrix for λ ∈ σ S, C . Then a generalized eigenvector x s corresponding to λ is a local minimum of 1.1 .
With grad x Φ x 2 Sx − λCx holds grad x Φ x s 0 and since the Hessian matrix S − λC of Φ is positive semidefinite, the vector x s is a minimal solution of Φ. Then, x s minimizes also F x subject to x t Cx d by the Lemma 3.8.
This value is corresponding in 3.11 with the inverse eigenvalue of problem 3.1 . Furthermore, it yields 3.18 So μ ≥ 0 follows, that is, λ * −1 I − R t −1 CR −1 is positive semidefinite, and for y ∈ R m we obtain By setting y Rx and with regular R, we get So, λ * λ E follows and x * s is a minimum of F x subject to x t Cx d.
Example 3.14. We minimize F : 1. So, we have d 1, S the identity matrix I 2 ∈ R 2×2 , and C ∈ R 2×2 a diagonal matrix with values −1 and 1. Then, we get the following generalized eigenvalue problem: with eigenvalues 1 and −1. Because of Theorem 3.1, we consider a generalized eigenvector α, 0 t with α ∈ R \ {0} only for λ 1. Then, 1, 0 t and −1, 0 t are solutions subject to ξ 2 1 − ξ 2 2 1. This result is conform to the geometric interpretation, since we are looking for x ξ 1 , ξ 2 t on the hyperbola ξ 2 1 − ξ 2 2 1 with minimal distance to the origin.

Reduction to an Eigenvalue Problem of Dimension rg C
In numerical applications, a generalized eigenvalue problem is mostly reduced to an eigenvalue problem, for example, by multiplication with S −1 . Thus, we obtain the inverse problem 3.14 from 3.1 see, e.g., 13 . But, S may be ill-conditioned, so that a solution of 3.14 may be numerical instable. Therefore, we present another reduction of 3.1 .
Many times, C is a sparse matrix with r : rank C ≤ rank S . This symmetric matrix C is diagonalizable in C P t DP with P orthogonal and D diagonal. Further, we assume that the first r diagonal entrees in D are different to 0. For the characteristic polynomial in 3.1 , it follows The order of p is r. We decompose these matrices in with S 1 , D 1 ∈ R r×r , S 2 ∈ R r× m−r , and S 3 ∈ R m−r × m−r . Now, we eliminate S 2 in PSP t by multiplications with Givens rotations G k ∈ R m×m , k 1, . . . , l, so that it follows: Journal of Applied Mathematics with Σ 1 , Δ 1 ∈ R r×r , Σ 2 , Δ 2 ∈ R m−r ×r , and Σ 3 ∈ R m−r × m−r . In 4.1 , we achieve with orthogonal G k , k 1, . . . , l

4.4
Because of p 0 det PSP t det S / 0 and p 0 det Σ 3 det Σ 1 , the submatrices Σ 1 , Σ 3 are regular and the generalized eigenvalues of det Σ 1 − λΔ 1 are different to zero. So, with y ∈ R r can be transformed to an equivalent eigenvalue problem with This system can be solved by finding the matrix X with Δ 1 X Σ 1 using the Gaussian elimination and determining the eigenvalues of X computing the QR algorithm 13 . Because all steps are equivalent, we have σ Δ −1 1 Σ 1 σ S, C , that is, the eigenvalues of 3.1 and 4.6 are the same.
With Theorem 3.13, we are looking for the smallest value λ * ∈ σ S, C and a corresponding generalized eigenvector x * s to minimize the problem 3.1 . So, yields. By substitution of y * s for Px * s , we obtain We decompose y * s into the subvectors y * s,r ∈ R r and y * s,m−r ∈ R m−r with y * s t y * s,r |y * s,m−r t .
Then, y * s,r is a generalized eigenvector for λ * of the problems 4.5 and 4.6 . Let y * s,r be an eigenvector to the smallest positive eigenvalue λ * of 4.6 . Since Σ 3 is regular, it follows in 4.8 and a generalized eigenvector x * s for λ * in 3.1 is given as:

4.10
Journal of Applied Mathematics 9

5.1
The equation f a, ξ, η 0 can also be written with x ξ, η t as The eigenvalues λ 1 , λ 2 of A characterize a conic uniquely 14 . Thus, we need 4λ 1 λ 2 4α 1 α 3 − α 2 2 > 0 for ellipses in f a, x 0. Furthermore, every scaled vector μa with μ ∈ R \ {0} describes the same zero-set of f. So, we can impose the constraint for ellipses with 4α 1 α 3 −α 2 2 1. For n n ≥ 6 given points ξ i , η i t ∈ R 2 , we want to find a parameter a ∈ R 6 , which This ellipse fitting problem is established and solved by Fitzgibbon 7 . With the following has exactly one positive solution λ * ∈ R. Because of Theorem 3.13 a corresponding generalized eigenvector a * to λ * minimizes the problem 5.5 and a * consists of the coefficients of an implicit given ellipse. A numerically stable noniterative algorithm to solve this optimization problem is presented by Halir and Flusser 15 . In comparison with Section 4, their method uses a special block decomposition of the matrices D and C.

Hyperbola Fitting
Instead of ellipses, O'Leary and Zsombor-Murray want to find a hyperbola to a set of scattered data x i ∈ R 2 9 . A hyperbola is a conic, which can uniquely be characterized by 4λ 1 λ 2 4α 1 α 3 − α 2 2 < 0 14 . So, we consider the constraint α 2 2 − 4α 1 α 3 1 and obtain the optimization problem:

Bookstein's Conic Fitting
In Bookstein's method, the conic constraint is restricted to where λ 1 , λ 2 are the eigenvalues of A in f 6 . There, it is λ 1,2 ∈ − √ 2/2, √ 2/2 and at least one of them different to 0. But the constraint 5.8 is not a restriction to a class of conics. Here, we determine an arbitrary conic, which minimizes The resulting data matrix D ∈ R n×6 is the same as for Fitzgibbon's problem. The constraint matrix C ∈ R 6×6 has a diagonal shape with the entrees 2, 1, 2, 0, 0, 0 , that is, all eigenvalues of C are nonnegative. In the case of a regular matrix S, the problem 5.9 is solved for a generalized eigenvector to the smallest value in σ S, C .

Approximation of Shrinked Aspheres
After the molding process in optical applications, the shrinkage of rotation-symmetric aspheres is implicitly defined for x ξ, ζ t in  subject to a t Ca 4r 2 ref .

5.13
This is also an application of 1.1 . The matrix C has the eigenvalues −2, 0, 1, and 2. So, the generalized eigenvalue problem in 3.1 with regular S D t D ∈ R 4×4 has two positive values in σ S, C . With Theorem 3.13, a generalized eigenvector a * ∈ R 4 to the smaller of both values solves 5.13 . The coefficients α i in the problems 5.5 and 5.13 correspond not to the same monomials ξ k ζ l . Hence, we have different matrices D and C.

Conclusion
In this paper, we present a minimization problem of least squares subject to absolute quadratic constraints. We develop a closed theory with the main result that a minimum is a solution of a generalized eigenvalue problem corresponding to the smallest positive eigenvalue. Further, we show a reduction to an eigensystem for numerical calculations. Finally, we study four applications about conic approximations. We analyze Fitzgibbon's method for direct ellipse-specific fitting, O'Leary's direct hyperbola approximation, Bookstein's conic fitting, and an optical application of shrinked aspheres. All these systems are attribute to the general optimization problem.