Loading [MathJax]/jax/output/HTML-CSS/jax.js

Friday, July 19, 2019

Quadratic Programming in Hilbert space. Part II. Adaptive algorithm.

This is a continuation to the previous post.

One disadvantage of the algorithm described in the previous post is that it requires an initial feasible point. Getting such initial point may require a significant amount of computations. Also in a case when number of the inequality constraints is large the rate of convergence may be slow. Often occurs a situation when an active set corresponding to the problem solution includes a small number of the inequality constraints, significantly less than their total number. In such cases it makes sense to apply a different approach for solving the problem (1),(4),(5). The essence of it is as follows: we construct a sequence of σk approximations, each of which is a solution to the auxiliary minimization problem of the functional (1) on a certain subset of constraints (4), (5). The auxiliary minimization problems are defined in such a way that the value of the functional (1) strictly increases for each subsequent approximation. These minimization problems are solved with the algorithm outlined in the previous post. Thus, in order to get solution of the original problem (1),(4),(5), it is necessary to solve a series of problems with a small number of constraints. Let's describe this in details.
Assume there is a subset of constraints (4), (5):
hi,φH=biiI1,hi,φHbiiIk2,Ik2I2 ,Ik=I1Ik2,
where I1 and I2 are defined in (7). Denote by Φk the convex set defined by constraints (28) and (29). Obviously Φ (a set defined by constraints (4) and (5)) is a subset of Φk.
Let ˆσk be a normal solution of system (28), (29). Assume that ˆσk is not satisfying all constraints (5), i.e. ˆσkΦ (otherwise ˆσk is the solution of the problem (1),(4),(5)), and itI2Ik2 is a number of the most violated constraint. We set
Ik+12=Ik2{it},Ik+1=I1Ik+12.
Let Φk+1 be the convex set defined by constraints with numbers from Ik+1, then Φk+1 is a subset of Φk:
Φk+1Φk.
and therefore the following inequality holds
ˆσk2H<ˆσk+12H.
This is the strict inequality because ˆσkΦk+1.
In a case when ˆσk+1 does not satisfy a constraint with a number from I2Ik+12 the described process is repeated. Once the inequality (30) always holds and the number of constraints is finite, the algorithm converges to the normal solution of the system (4),(5) in a finite number of steps.
The computational details of the algorithm are described below.
Let ˆAk be an active set corresponding to the normal solution ˆσk of system (28), (29):
ˆAk={iIk:hi,ˆσkH=bi},
then
ˆσk=jˆAkμjhj,jˆAkgijμj=bi ,iˆAk,gij=hi,hjH ,iˆAk,  jˆAk.
For every iI2Ik2 we compute values
ˆeki=hi,ˆσkHbi=jˆAkgijμjbi,
 and find among them the largest, whose index is denoted as it, i.e.:
ˆekit=maxiI2Ik2{ˆeki}.
If  ˆekit0, then ˆσk satisfies to all constraints (5) and it is a solution of the original problem (1), (4), (5), so σ=ˆσk. Otherwise a new auxiliary problem, which consists in minimizing the functional (1) on solutions of the augmented constraint subsystem 
hi,φH=biiI1,hi,φHbiiIk2,hit,φHbit,Ik+12=Ik2{it} ,Ik+1=I1Ik+12,
should be solved. An initial feasible point required to solve this problem can be defined as the normal solution ˜σk+1 of the system
hi,φH=biiˆAk,hit,φH=bit,
Let ˜Ik+1=ˆAk{it} then ˜σk+1 is presented as 
˜σk+1=j˜Ik+1μjhj ,j˜Ik+1gijμj=bi ,i˜Ik+1 .
Since matrix of the system (34) differs from the matrix of the system (32) in only the last row and column, its Cholesky decomposition can be obtained by modifying the Cholesky decomposition of the system (32) matrix as it is described in Algorithms for updating a Cholesky factorization post.
Initial set I0 can be defined as follows:
I0=I1{i0,i0+M},
where i0 is such that
δi0=minL+1iSδi.

The idea of this algorithm was originally published in the paper [2] where it was applied for constructing of "splines in convex set". In that case, the minimizing objective function of the related variational problem was a semi-norm, and existing and uniqueness of the problem solution was not trivial. Corresponding auxiliary quadratic programming problems were not strictly convex ones and for their numerical solution a variant of coordinate descent algorithm was used. The algorithm showed good results when solving one-dimensional problems.

Some numerical results and comparing of the described algorithms will be presented in the next post.

References

[1] I. Kohanovsky, Data approximation using multidimensional normal splines. Unpublished manuscript. 1996.
[2] A. Kovalkov, On an algorithm for the construction on splines with discrete constraints of the type of inequalities, Serdica Bulgariacae mathematicae publicationes, Vol.9, N.4, 1983