Méthode du gradient conjugué

Gentil · October 2019

Bonjour à la communauté
J'ai essayé d'écrire un code pour le gradient conjugué tel que je l'ai compris mais il ne marche pas vraiment. Je rappelle qu'il s'agît de trouver le minimum de $$
\frac{1}{2} \langle Ax,x \rangle + \langle b,x \rangle

$$ Ici une subtilité est que $A$ est une fonction linéaire et non une matrice. Mais l'idée reste la même.
Je pars de $x_{n} = 0$. Je prends la direction opposée du gradient qui est $$
d_{n+1} = -(Ax_{n}-b)
$$ puis je transforme cette direction $d_{n+1}$ en un vecteur unitaire orthogonal pour $A$ par Gram Smith. Et suivant ce vecteur je suis censé arriver au minimum.

def solveGC(A, b):
    d = len(b)
    e = np.zeros( (d+1,d) )
    dir = 0 
    sol = np.zeros(d)
    for n in range(d+1):
        dir = -A(sol) + b
        
        proj = 0
        if n >0:
            for i in range(n):
                proj += np.dot( dir.T, A(e[ i].T) )*e[ i].T
        
        if la.norm( np.dot( (dir-proj).T, A(dir-proj) )  ) >0:
            e[n] = (dir-proj)/np.dot( (dir-proj).T, A(dir-proj) ) 
        else:
            e[n] = 0
        sol = A(e[n].T)
    return sol

Gentil · October 2019

J'ai modifié un peu mon code car je ne tenais pas compte de certains points.

[Mettre une ' ' entre [ et i sinon, l'afficheur du forum le prend pour une bannière BBcode de mise en italique. AD]

def solveGC(A, b, tol=1e-6):
    d = len(b)
    # initialisation base orthogonale des directions
    e = np.zeros( (d+1,d) )
    # Point de départ : 0
    xk = np.zeros(d)
    # initialisation solution
    sol = np.zeros(d)
    # initialisation direction opposée au gradient
    dk = np.zeros(d)

    for n in range(1,d+1):
        #Direction opposée au gradient
        dk = b-A(xk)

        # On calcul le minimum de F dans la direction opposée au gradient : F(xk + t*dk)
        min = - np.dot(A(xk).T,dk)/np.dot(A(dk).T,dk) # On divise par 0 ? La direction dans le noyau ?
        # Non car direction famille libre pour Gram Schmidt

        #Calcul du prochain point
        xk = xk + min*dk
        # Calcul de la prochaine direction
        proj = 0

        if n >0:
            for i in range(n):
                proj += np.dot( dk.T, A(e[ i].T) )*e[ i].T
        else:
            e[0] = dk/np.dot(A(dk).T,dk)
        e[n] = (dk-proj)/np.dot( (dk-proj).T, A(dk-proj) )
        sol = A(e[n].T)
    return sol

Méthode du gradient conjugué

Réponses

Bonjour!

Catégories

In this Discussion

Qui est en ligne 7