Ellipsoïde en boite

En cherchant à exprimer des intervalles de confiance pour des coefficients lors d’une régression linéaire multiple, j’ai approché ma statistique de khi-deux (χ²) par une fonction quadratique convexe au voisinage de son minimum. L’ensemble des vecteurs de paramètres en dessous du seuil de khi-deux acceptable est alors un ellipsoïde (ou plus exactement sa généralisation en dimension quelconque) dont je cherche la projection sur chaque axe. Les sections par les sous-espaces engendrés par des axes de coordonnées se traitent par des formules analogues. Ces formules permettent de construire ensuite des algorithmes d’optimisation (de la fonction ou des projections de l’ellipsoïde) sous la contrainte que les vecteurs soient de coordonnées toutes positives.

Projections sur les axes

On considère un ellipsoïde (plein) de centre C pouvant être décrit par une inéquation de la forme (X − C)⊤S(X − C) ≤ z où S est une matrice symétrique définie positive de taille n et z est un réel positif.

Ellipse et gradient aux points optimaux pour la projection sur l’axe horizontal. Ces points optimaux ne sont pas les sommets de l’ellipse.

Sa projection sur un axe de coordonnée est un segment dont les extrémités sont obtenues en optimisant la coordonnée Xi sur le bord de l’ellipsoïde. La méthode du multiplicateur de Lagrange montre qu’en un tel point, le gradient de la fonction quadratique f(X) = (X − C)⊤S(X − C) doit être colinéaire au vecteur élémentaire Ei.

Comme ce gradient s’écrit ∇f = 2S(X − C), une condition nécessaire de l’optimisation est l’existence d’un réel λ tel que X = C + λS−1Ei. Un tel point sera dit optimal pour la coordonnée Xi.

Dans ce cas, on trouve f(X) = λ2(S−1)i,i d’où λ = ±√(z/(S−1)i,i).

Chaque coordonnée du vecteur X est donc encadrée par les valeurs Xi± = Ci ± √(z(S−1)i,i).

Restriction à un sous-espace

Soit A ⊂ ⟦1, n⟧. On note RA le sous-espace engendré par les axes indexés par les éléments de A. Si on cherche à encadrer les composantes d’une restriction de l’ellipsoïde à ce sous-espace  on peut introduire la matrice J représentant l’injection canonique de RA dans Rn. On réécrit alors X = JX′ puis on pose S′ = J⊤SJ (matrice symétrique définie positive représentant la restriction de la forme quadratique à F) et C′ = S′−1J⊤SC pour obtenir :

(X′ − C′)⊤S′(X′ − C′) = z − C⊤SC + C′⊤S′C′

En notant z′ le second membre de cette égalité, on retrouve les expressions Xi′± = C′i ± √(z′(S′−1)i,i), à condition bien sûr que z′ soit positif, ce qui n’est pas forcément le cas : si z′ < 0, l’ellipsoïde est disjoint du sous-espace RA.

En particulier, si A = {i}, on a S′ = Si,i donc C′ = 1/Si,iSiC

Point de coordonnées positives

On s’intéresse maintenant uniquement à la partie de l’ellipsoïde dont les points ont des coordonnées toutes positives, c’est-à-dire son intersection avec l’orthant positif. La première question qui se pose est l’existence d’un point dans cette intersection, et sa détermination.

L’algorithme qui suit est une méthode d’ensemble actif (active set) décrite par exemple par J. Nocedal et S. W. Wright dans Numerical Optimization, §16.5 Active-Set Methods for Convex QPs, Springer 2006. Il est composé d’une alternance de deux phases, que j’appelle restriction et extension, définissant une suite de points de valeur strictement croissante selon f, chacun étant situé au centre d’une restriction de l’ellipsoïde à un sous-espace engendré par des axes de coordonnées. L’ensemble de ces sous-espaces étant fini (indexé par l’ensemble des parties de ⟦1, n⟧), l’algorithme termine nécessairement.

L’initialisation consiste à tester la positivité des coordonnées du centre (ce qui achèverait immédiatement la recherche), sinon on définit un point M0 par projection de ce centre sur l’orthant positif, autrement dit en annulant les coordonnées négatives du centre, et on définit le sous-espace F0 = Rn.

Phase de restriction en deux étapes. Le point Mk+1 est bloqué sur le segment entre Mk et Ck par l’annulation de la deuxième coordonnée. Le point Cn+1 est le milieu du segment intersection de l’ellipse pleine avec l’axe des abscisses. La droite (CkCk+1) n’est pas l’axe de l’ellipse.

La phase de restriction est une boucle associée à une suite strictement décroissante (pour l’inclusion) de sous-espaces Fk, ce qui assure là encore sa terminaison. À chaque étape, on teste la positivité des coordonnées du centre Ck de la restriction de l’ellipsoïde à Fk. Le cas échéant, la boucle se termine, sinon, on relie le point Mk à ce centre et le nouveau point Mk+1 est le plus avancé sur ce segment qui reste à coordonnées positives. Le nouveau sous-espace Fk+1 est obtenu en annulant les coordonnées qui bloquent l’avancée le long du segment.

Par stricte convexité de la fonction f, ces valeurs décroissent strictement le long de chaque segment [MkCk] en direction du centre et par conséquent le long de la suite des points (Mk).

La phase d’extension consiste à étendre le sous-espace selon une coordonnée négative du gradient. Le centre de la restriction de l’ellipsoïde sur ce nouveau sous-espace est alors nécessairement avec une coordonnée positive dans cette direction, et on peut relancer une phase de restriction.

Lorsque toutes les composantes du gradient sont positives, cela signifie que le point minimise la fonction f sur l’orthant positif, ce qui arrive en temps fini. À chaque étape de calcul, on peut arrêter l’algorithme si la fonction est évaluée en dessous du seuil de définition de l’ellipsoïde. Si à la fin de l’algorithme, la minimisation de la fonction reste au dessus du seuil, c’est que l’ellipsoïde est disjoint de l’orthant positif.

Projections de la partie positive

Recherche des points optimaux pour la projection sur l’axe des abscisses. Ici deux étapes terminent la phase de restriction pour la maximisation, tandis que la minimisation aboutit à une borne nulle.

Pour évaluer les projections sur les axes de l’intersection de l’ellipsoïde avec l’orthant positif, on peut appliquer un algorithme analogue au précédent, mais cette fois le point de départ est un point de coordonnées positives dans l’ellipse, et ses successeurs sont construits en direction des points optimaux au lieu de viser les centres.

Par unicité des points optimaux (pour chaque direction d’axe de coordonnée et chaque sens), la coordonnée à optimiser est strictement améliorée à chaque étape d’une phase de restriction. Dans le cas d’une minimisation, il se peut qu’un des points d’étape ait une coordonnée nulle, auquel cas la recherche s’arrête et 0 est bien la valeur minimale à trouver.

Lors de la phase d’extension, une coordonnée à maximiser est toujours associée à une composante positive du gradient, mais une coordonnée à minimiser est associée à une composante négative du gradient (si la minimisation est elle-même différente de 0). On étend l’espace de recherche dans la direction d’une coordonnée négative du gradient à l’exception de la coordonnée à minimiser.

Programmation en Python

Le programme ci-dessus définit une classe Ellipsoide avec les méthodes suivantes :

projections
Bornes des coordonnées des points de l’ellipsoïde
restriction
Crée ou récupère les éléments de la restriction mise en cache
relaxe
Complète les coordonnées manquantes d’un vecteur restreint par des zéros
point_positif
Détermine un point de l’ellipsoïde avec des coordonnées positives en utilisant une méthode d’ensemble actif (active-set)
proj_positives
Bornes des coordonnées des points de l’intersection de l’ellipsoïde avec l’orthant positif

Une même restriction pouvant servir plusieurs fois dans la fonction proj_positives, ses éléments sont enregistrés en un objet de la sous-classe Restriction dans une table indexée par des chaines de 0 et 1 (codant le masque de l’ensemble actif).

#!/usr/bin/python3
# -*- coding: utf-8 -*-

	
import numpy
	
class Ellipsoide(object):
	"""Définit un ellipsoïde plein de centre C, de matrice S et de niveau z
	dont les points M vérifient l’inégalité (M-C)^T S (M-C) ≤ z"""
	def __init__(self, centre, matrice, niveau):
		assert niveau >= 0
		self.C = numpy.array(centre)
		self.S = numpy.array(matrice)
		self.Sinv = numpy.linalg.inv(matrice)
		self.z = niveau
		self.restrictions = {}
	def projections(self):
		"""Bornes des coordonnées des points de l’ellipsoïde"""
		return [(c-r, c+r) for (c, r) in ((self.C[i], (self.z*self.Sinv[i,i])**0.5) 
			for i in range(self.C.shape[0]))]
	class Restriction(object):
		"""Définit les éléments de la restriction de la forme quadratique à l’espace
		engendré par les axes correspondant aux éléments de A"""
		def __init__(self, parent, A):
			if A:
				S = numpy.array([[parent.S[j,k] for j in A] for k in A])
				self.Sinv = numpy.linalg.inv(S)
				H = numpy.array([parent.SC[k] for k in A])
				self.C = numpy.linalg.solve(S, H)
				self.z = parent.z0 + H @ (self.Sinv @ H)
				self.indices = {k:j for (j, k) in enumerate(A)}
			else:
				self.C = numpy.array(0)
				self.z = parent.z0
		def optimal(self, i, epsilon):
			"""Point optimal pour la coordonnée i dans le sens du signe d’epsilon"""
			j=self.indices[i]
			return self.C + epsilon*(self.z/self.Sinv[j,j])**0.5*self.Sinv[j,:]
	
	def restriction(self, A):
		"""Crée ou récupère les éléments de la restriction mise en cache"""
		grille = ["0" for _ in range(self.C.shape[0])]
		for j in A:
			grille[j] = "1"
		key = "".join(grille)
		if not self.restrictions:
			self.SC = self.S @ self.C
			self.z0 = self.z-self.C @ (self.SC)
		if key not in self.restrictions:
			self.restrictions[key] = Ellipsoide.Restriction(self, A)
		return self.restrictions[key]
	def relaxe(X, A):
		"""Complète les coordonnées manquantes (hors de A) par des zéros"""
		Y = numpy.zeros(self.C.shape[0])
		for (j,k) in enumerate(A):
			Y[k] = X[j]
		return Y
	def point_positif(self):
		"""Détermine un point de l’ellipsoïde avec des coordonnées positives
		en utilisant une méthode d’ensemble actif (active-set)"""
		if numpy.amin(self.C)>=0:
			return self.C
		M = numpy.maximum(self.C, 0)
		V = M-self.C
		y = V @ (self.S @ V)
		if y <= self.z:
			return M
		n = self.C.shape[0]
		A = [i for i in range(n) if self.C[i]>=0]
		M = numpy.array([M[j] for j in A])
		while 1: 
			# phase de restriction
			r = self.restriction(A)
			while numpy.ndim(r.C) and numpy.amin(r.C)<0:
				t = 1
				for j in range(M.shape[0]):
					if r.C[j]<0 and (1-t)*M[j]+t*r.C[j]<0:
						t = M[j]/(M[j]-C[j])
						jmin = j
				M = numpy.maximum((1-t)*M+t*r.C, 0)
				M[jmin] = 0
				V = M-r.C
				y = V @ (r.S @ V)
				if y <= r.z:
					return self.relaxe(M, A)
				A.pop(jmin)
				M = numpy.delete(M, jmin)
				r = self.restriction(A)
			if r.z >= 0:
				return self.relaxe(r.C, A)
			# phase d’extension
			G = sum(x * self.S[k,:] for (x,k) in zip(r.C, A)) - self.SC
			j = 0
			seuil = 0
			i0 = j0 = -1
			for i in range(n):
				if j<len(A) and A[j]==i:
					j=j+1
				elif G[i] < seuil:
					i0, j0 = i, j
					seuil = G[i]
			if i0<0:
				raise NameError("Pas de point de coordonnées positives")
			A.insert(j0,i0)
			M = numpy.insert(M, j0, 0)
	
	def proj_positives(self, M0=None):
		"""Bornes des coordonnées des points de l’intersection de l’ellipsoïde
		avec l’orthant positif"""
		resultat = []
		n=self.C.shape[0]
		for i in range(n):
			bornes = []
			V = (self.z/self.Sinv[i,i])**0.5*self.Sinv[i,:]
			for epsilon in (-1, 1):
				X = self.C + epsilon*V
				if numpy.amin(X)>=0:
					bornes.append(X[i])
				else:
					if M0 is None:
						M0 = self.point_positif()
					M = M0
					A = list(range(n))
					recherche = True
					while recherche:
						# phase de restriction
						while numpy.amin(X)<0:
							t = 1
							for j in range(M.shape[0]):
								if X[j] < 0 and (1-t)*M[j]+t*X[j] < 0:
									t = M[j]/(M[j]-X[j])
									jmin = j
							if A[jmin] == i and epsilon==-1:
								recherche = False
								bornes.append(0)
								break
							M = numpy.maximum((1-t)*M+t*X, 0)
							A.pop(jmin)
							M = numpy.delete(M, jmin)
							X = self.restriction(A).optimal(i, epsilon)
						# phase d’extension
						if recherche:
							G=sum(x*self.S[k,:] for (x,k) in zip(X, A))-self.SC
							j=0
							for k in range(n):
								if j<len(A) and A[j]==k:
									j=j+1
								elif G[k]<0:
									A.insert(j,k)
									M = numpy.insert(M, j, 0)
									X = self.restriction(A).optimal(i, epsilon)
									break
							else:
								recherche = False
								bornes.append(X[A.index(i)])
			resultat.append((bornes[0], bornes[1]))
		return resultat