Discussion :

[Alexis.M]

Bonjour,

Je cherche à résoudre efficacement le problème suivant:

(a * I_n + QQ^T) x = y

Où a est un scalaire positif, I_n est la matrice identité de taille n x n, et Q une matrice n x m telle que m << n. x et y sont des vecteurs de taille n.

L'inverse de la matrice est a priori calculable par:

(a * I_n + QQ^T)^(-1) = 1/a * ( I_n + Q (a * I_m + Q^TQ)^(-1)Q^T))

Ce qui permet de beaucoup réduire les temps de calcul puisque la matrice à inverser est de taille m x m au lieu de n x n. Cependant, j'imagine que, comme pour la solution de Ax = y, le calcul de l'inverse n'est pas la méthode la plus adaptée pour des problèmes de précision numérique et de temps de calcul.

Quelle est donc la meilleure façon de résoudre ce type de système ?

D'avance merci,

Alexis

[bertry]

Salut,

le calcul de l'inverse n'est pas la méthode la plus adaptée pour des problèmes de précision numérique et de temps de calcul.

Habituellement, une décomposition LU ou une décomposition de Cholesky est utilisée pour résoudre les systêmes linéaires.

[Alexis.M]

Habituellement, une décomposition LU ou une décomposition de Cholesky est utilisée pour résoudre les systêmes linéaires.

Oui mais mon cas a une structure très particulière. Il n'y a pas d'intérêt à calculer naïvement la décomposition de Cholesky de (a*I_n + QQ^T), à moins de disposer d'une méthode spécialisée pour la structure du problème, faute de quoi on se retrouve avec la complexité O(n^3) habituelle. Alors que la formule de l'inverse montre qu'il doit être possible de calculer le tout en O(m^3 + n^2m). Avec m << n, le gain en temps de calcul est considérable.

[souviron34]

Si je comprend bien, tu as un système sur-dimensionné, c'est ça ?

N équations à m inconnues ?

[acx01b]

A x = b

où A a plus de colonnes que de lignes (donc x plus grand que b)

supposons que A est de rang plein (pas de lignes dépendantes)
x = x1 + x2
où x1 = A^T y est une combinaison linéaire des lignes de A,
et A x2 = 0 (donc x2 appartient à ker(A))

tu résous pour y (A A^T est inversible) :
A A^T y = b

puis si tu veux tu cherches une base de ker(A) pour trouver l'ensemble des x2 possibles

si A n'est pas de rang plein c'est pareil mais A A^T n'est pas inversible donc il peut ne pas y avoir de solution pour y

pour trouver ker(A) je pense à Gram-Schmidt

concernant la régularisation, c'est A A^T y = b qu'il faut régulariser

[souviron34]

ben moi je pensais au Simplex, si effectivement il ya plus d"équations que d'inconnues

[acx01b]

ben moi je pensais au Simplex, si effectivement il ya plus d"équations que d'inconnues

quel rapport avec le simplexe ? il n'y a pas de contrainte (linéaire) ici

le simplexe c'est dans le cas où il y a des contraintes linéaires ce qui assure que la solution du truc linéaire à minimiser est sur un des sommets du polytope acceptable (on oublie les cas dégénérés mais bon l'idée c'est ça)

[acx01b]

A x = b

si maintenant x est plus petit que b,
c'est la même idée

tu décomposes b = b1+b2 où A^T b2 = 0
si b2 n'est pas égal à 0 alors il n'y a pas de solution

sinon

A^T A x = A^t b

que tu résous pour x (et où tu peux régulariser A^T A)

[souviron34]

quel rapport avec le simplexe ? il n'y a pas de contrainte (linéaire) ici

Je posais la question - je ne comprend pas exactement la formulation. Un Simplex et applicable sur un systèle linéaire, dans le cas où on a plus d'équations (n) que d'inconnues (m) (ce qui revient à avoir des contraintes en nombre n-m)..

[acx01b]

hors sujet sur le simplexe

regarde wiki/Algorithme_du_simplexe

la forme la plus couramment utilisée pour présenter et étudier les algorithmes, qui est celle supposée par l'algorithme du simplexe révisé, est la forme standard :

min c^T x
Ax = b
et x >= 0

où c,b sont des vecteurs connus, A une matrice connue, et x un vecteur qu'on cherche

la principale difficulté est de montrer que x* (le x optimal) est forcément sur un sommet du polytope admissible

à mon avis le plus simple c'est d'abord de montrer que soit x = infini sur une des composante soit le peut être vu comme borné (en rajoutant des contraintes u^T x <= v qui ne changent rien à la solution), puis de dire que si x* est dans l'intérieur du polytope, alors x* est sur le segment [ab] où a est un sommet et b un point d'une face, a et b sont admissibles et c^Tx est une fonction linéaire donc
c^Ta < c^Tb ou c^Tb < c^Ta ou c^Ta = c^Tb
dans les deux premiers cas x* n'est pas la solution et dans le troisième a et b sont également solution

si x* est sur une face c'est le même principe


s'en suit l'algorithme du simple où l'on se balade de sommets en sommets adjacents

[Alexis.M]

Si je comprend bien, tu as un système sur-dimensionné, c'est ça ?

N équations à m inconnues ?

Non pas du tout, il s'agit d'un système linéaire avec autant d'équations que d'inconnues.

La matrice A = (a * I_n + QQ^T) est de taille n x n, symétrique et de rang plein, mais la composante Q de taille n x m, est elle de faible rang.

Ceci a pour conséquence que, lorsque je calcule A^(-1), je peux utiliser le théorème de l'inverse comme dans mon premier message pour réduire les calculs. Il est donc presque certain qu'il doit y a voir moyen de résoudre Ax = b de manière plus efficace et plus précise qu'en passant par le calcul de l'inverse.

[Alexis.M]

Je pense avoir trouvé une solution raisonnable. Puisque

x = (a * I_n + QQ^T)^(-1) y
x = 1/a * (I_n + Q(a * I_m + Q^TQ)^(-1)Q^T) y # Formule de Woodbury
x = 1/a * (y + Q(a * I_m + Q^TQ)^(-1)Q^T y )

Je peux résoudre:
(a*I_m + Q^TQ) z = Q^Ty
puis:
x = 1/a * (y + Qz)

Cela permet notamment d'utiliser la décomposition de Cholesky pour calculer z. Il faut:

• (m * N) opérations pour calculer Q^Ty
• (m^2 * N) opérations pour Q^TQ
• O(m^3) pour calculer z
• (N + N * m) pour calculer (y + Qz)

Soit une complexité asymptotique à la louche en O(m^3 + m^2 * N) à comparer aux O(N^3) de la solution naïve.