opti-non-lin/1-coordinate-descent.ipynb

In [1]:

import numpy as np

In [2]:

# Un premier exemple de descente par coordonnées

# Itéré intial
x = np.array([-1/2, -1])

maxiter = 10
for k in range(maxiter):
    #On mets (de façon optimale) x1
    x[0] = (6*x[1]+4)/10
    #On mets (de façon optimale) x2
    x[1] = (6*x[0]-3)/10
    print(f"A l'itété {k}, on trouve x={x}")

A l'itété 0, on trouve x=[-0.2  -0.42]
A l'itété 1, on trouve x=[ 0.148  -0.2112]
A l'itété 2, on trouve x=[ 0.27328  -0.136032]
A l'itété 3, on trouve x=[ 0.3183808  -0.10897152]
A l'itété 4, on trouve x=[ 0.33461709 -0.09922975]
A l'itété 5, on trouve x=[ 0.34046215 -0.09572271]
A l'itété 6, on trouve x=[ 0.34256637 -0.09446018]
A l'itété 7, on trouve x=[ 0.34332389 -0.09400566]
A l'itété 8, on trouve x=[ 0.3435966  -0.09384204]
A l'itété 9, on trouve x=[ 0.34369478 -0.09378313]

In [3]:

def CDquadratic(Q,c,p,x,maxiter=10**2):
    # Ce code applique la méthode de descente par coordonnées
    # au problème min_x 0.5 x'*Q*x-c'*x+p
    # Le x fournit en paramètre est l'itéré initial
    # Q nous est donnée sous forme symétrique
    n = len(x)
    # Une itération coûte O(n^2)
    for k in range(maxiter):
        # Dans cette boucle, on mets les variables à jour
        for i in range(n):
            somme = 0
            for j in range(n):
                if not j==i:
                    somme += Q[i,j]*x[j]
            x[i] = (c[i]-somme)/Q[i,i]
        print(f"A l'itété {k}, on trouve x={x}")
    return x

In [4]:

Q = np.array([[10, -6],[-6, 10]])
c = np.array([4, -3])
p = 0
x = np.array([-1/2, -1])
xopt = CDquadratic(Q,c,p,x,10)

A l'itété 0, on trouve x=[-0.2  -0.42]
A l'itété 1, on trouve x=[ 0.148  -0.2112]
A l'itété 2, on trouve x=[ 0.27328  -0.136032]
A l'itété 3, on trouve x=[ 0.3183808  -0.10897152]
A l'itété 4, on trouve x=[ 0.33461709 -0.09922975]
A l'itété 5, on trouve x=[ 0.34046215 -0.09572271]
A l'itété 6, on trouve x=[ 0.34256637 -0.09446018]
A l'itété 7, on trouve x=[ 0.34332389 -0.09400566]
A l'itété 8, on trouve x=[ 0.3435966  -0.09384204]
A l'itété 9, on trouve x=[ 0.34369478 -0.09378313]

In [5]:

def leastsquares(A,b,x,maxiter=10**2):
    Q = 2*A.T@A
    c = 2*A.T@b
    p = b.T@b
    x = CDquadratic(Q,c,p,x,maxiter)
    return x

In [6]:

A = np.array([[1,0], [0,1], [1,1]])
b = np.array([3,2,-3])
x = np.array([1.0, 0.0]) # Attention au entiers
xopt = leastsquares(A,b,x,10)

A l'itété 0, on trouve x=[ 0.  -0.5]
A l'itété 1, on trouve x=[ 0.25  -0.625]
A l'itété 2, on trouve x=[ 0.3125  -0.65625]
A l'itété 3, on trouve x=[ 0.328125  -0.6640625]
A l'itété 4, on trouve x=[ 0.33203125 -0.66601562]
A l'itété 5, on trouve x=[ 0.33300781 -0.66650391]
A l'itété 6, on trouve x=[ 0.33325195 -0.66662598]
A l'itété 7, on trouve x=[ 0.33331299 -0.66665649]
A l'itété 8, on trouve x=[ 0.33332825 -0.66666412]
A l'itété 9, on trouve x=[ 0.33333206 -0.66666603]