Otimização/Método de gradientes conjugados

Predefinição:Wikipedia

Algumas considerações históricas

• Este método foi originalmente proposto por Hestenes e Stiefel, em 1952.
• Seu objetivo inicial foi a resolução de problemas quadráticos sem restrições, mas logo o mesmo foi estendido para casos mais gerais.

O método

Este método pode ser considerado sob dois pontos de vista:

• Como um método de descida, com busca linear exata;
• Como um método de resolução de sistema linear, baseado em um processo de ortogonalização.

Exemplos de conjuntos convexos e côncavos

• 
  Este é um conjunto convexo, pois todo segmento com extremidades no conjunto está totalmente contido no conjunto.
• 
  Este é um conjunto côncavo, pois existe um segmento com extremidades no conjunto que não está totalmente contido no conjunto.

Predefinição:Definição

Nota: Uma matriz é definida positiva se, e somente se, todos os seus autovalores são positivos.

Tem-se:

$\mathrm{\nabla }f:{ℝ}^n\mapsto ℝ$

${\mathrm{\nabla }}^{2}f:{ℝ}^n\mapsto ({ℝ}^n{)}^2$

Sendo $f(x)=\frac{1}{2}{x}^{\mathrm{\top }}Ax+a^{\mathrm{\top }}x+\alpha$, segue em particular que $\mathrm{\nabla }f=Ax+a$ e ${\mathrm{\nabla }}^{2}f=A=P^{\mathrm{\top }}\mathrm{\Lambda }P$, onde $\mathrm{\Lambda }$ é uma matriz diagonal cujos elementos da diagonal são os autovalores de $A$ e $P$ é uma matriz onde as colunas são os autovetores correspondentes aos autovalores.

Note que $A$ é uma matriz simétrica, pois é a matriz Hessiana de uma função com segundas derivadas parciais contínuas, e consequentemente vale $\frac{{\partial }^{2}f}{\partial x\partial y}=\frac{{\partial }^{2}f}{\partial y\partial x}$.

Para introduzir o método de direções conjugadas, serão consideradas somente funções quadráticas.

Uma condição necessária de primeira ordem para que $x$ seja um ponto de mínimo para a função $f$ é que $\mathrm{\nabla }f(x)=0$. Para o presente caso, a função $f$ é convexa, então, a condição necessária $\mathrm{\nabla }f(x)=0$ também é suficiente.

No caso de uma função quadrática, tem-se $\mathrm{\nabla }f(x)=0\iff Ax+a=0$, ou seja, $x$ é solução do sistema linear $Ax=-a$.

A resolução de um sistema linear nem sempre pode ser feita numericamente de forma eficiente. Por exemplo, se a matriz do sistema é:

$A=\left[\begin{array}{cc}10^{-20}& 1\\ 1& 10^{20}+1\end{array}\right]$

A solução do sistema linear corresponde à interseção entre duas retas quase paralelas, e os erros de truncamento podem causar imprecisão na solução obtida computacionalmente.

Analiticamente, o sistema $Ax=-a$ tem $x=-A^{-1}a$ como solução. Então alguém poderia se perguntar: qual o problema em resolver esse sistema linear, se basta calcular a inversa da matriz $A$ e multiplicar pelo vetor $-a$? A resposta é que o calculo da inversa de uma matriz em geral é impraticável computacionalmente, por ter custo muito alto. Por isso, nas situações práticas, onde as matrizes tem ordem bem maior do que 2 (digamos 1000), o cálculo de matrizes inversas não é uma opção.

Assim, com o intuito de desenvolver um método computacional para o cálculo de minimizadores, é preciso utilizar outras técnicas. Considere o seguinte:

Em um método de descida tem-se sempre uma sequencia $\{x_k,t_k,d_k\}\in \mathbb{N}$, com $x_{k+1}=x_k+t_k{d}_k$ e $t_k$ é um minimizador de $f(x_k+t{d}_k):t\in ℝ$

$\mathrm{\nabla }f(x)=Ax+a$

e

$0=\mathrm{\nabla }f(x_{k+1})=A{x}_{k+1}+a=A(x_k+t_k{d}_k)+a=A{x}_k+t_{k}A{d}_k+a$

Logo, $t_{k}A{d}_k=-(A{x}_k+a)$ e multiplicando por $d_k^{\mathrm{\top }}$ obtem-se $t_k{d}_k^{\mathrm{\top }}A{d}_k=-(d_k^{\mathrm{\top }}A{x}_k+d_k^{\mathrm{\top }}a)$. Consequentemente, o valor de $t_k$ é dado por

$t_k=\frac{-(d_k^{\mathrm{\top }}A{x}_k+a^{\mathrm{\top }}{d}_k)}{d_k^{\mathrm{\top }}A{d}_k}$

Deste modo, o método consistirá de escolher em cada etapa $k$ uma direção $d_k$, e calcular o coeficiente $t_k$ pela fórmula anterior, para gerar o próximo ponto $x_{k+1}$. Mas como escolher a direção $d_k$?

Dado $x_k$ e escolhido $d_k$, defina $\theta :R\mapsto ℝ$ como $\theta (t)=f(x_k+t{d}_k)$, ou seja, $\theta$ é a restrição da função $f$ à reta que passa pelo ponto $x_k$ e que tem direção $d_k$. Logo, derivando a expressão de $\theta$ em relação a $t$, obtem-se

${\theta }^{\prime }(t)=\mathrm{\nabla }f(x_k+t{d}_k{)}^{\mathrm{\top }}{d}_k$

Então, no ponto de mínimo, $x_{k+1}$, tem-se

$0=\mathrm{\nabla }f(x_{k+1}{)}^{\mathrm{\top }}{d}_k$

Ou seja, a direção $d_k$ a ser seguida a partir do ponto $x_k$ é ortogonal ao gradiente da função $f$, no ponto $x_{k+1}$.

Esquema do método de descida

$x_{k+1}=x_k+t_k{d}_k=(x_{k-1}+t_{k-1}{d}_{k-1})+t_k{d}_k=\dots =x_1+t_1{d}_1+\dots +t_k{d}_k=x_1+\sum _{i=1}^k{t}_i{d}_i$

Seja $\overline{x}$ o minimizador da função $f$. Tem-se

$x_{k+1}-\overline{x}=x_1-\overline{x}+\sum _{i=1}^k{t}_i{d}_i$

Mas $0=\mathrm{\nabla }f(\overline{x})=A\overline{x}+a$ implica que $a=-A\overline{x}$, logo

$\mathrm{\nabla }f(x)=Ax+a=Ax-A\overline{x}=A(x-\overline{x})$

e consequentemente

$A(x_{k+1}-\overline{x})=A(x_1-\overline{x})+\sum _{i=1}^k{t}_{i}A{d}_i$

Donde $\mathrm{\nabla }f(x_{k+1})=\mathrm{\nabla }f(x_1)+\sum _{i=1}^k{t}_{i}A{d}_i$. Portanto $0=\mathrm{\nabla }f(x_{k+1}{)}^{\mathrm{\top }}{d}_k=\mathrm{\nabla }f(x_1{)}^{\mathrm{\top }}{d}_k+\sum _{i=1}^k{t}_i{d}_i^{\mathrm{\top }}A{d}_k$.

Usando o resultado desse exercício, tem-se ainda que $0=\mathrm{\nabla }f(x_{k+1}{)}^{\mathrm{\top }}{d}_k=\mathrm{\nabla }f(x_1{)}^{\mathrm{\top }}{d}_k+\sum _{i=1}^k{t}_i(B{d}_i{)}^{\mathrm{\top }}(B{d}_k)$

Fazendo $\delta =Bd$, o método do gradiente conjugado escolhe as direções de descida tais que ${\delta }_i^{\mathrm{\top }}{d}_j=0,\mathrm{\forall }i≠j$. Mas quando $i≠j$, tem-se na expressão apresentada anteriormente apenas $0=\mathrm{\nabla }f(x_1{)}^{\mathrm{\top }}{d}_k+t_k(B{d}_k{)}^{\mathrm{\top }}(B{d}_k)=\mathrm{\nabla }f(x_1{)}^{\mathrm{\top }}{d}_k+t_k{d}_{k}A{d}_k$

Finalmente, tem-se o algoritmo para este método.

Algoritmo de Hestenes-Stiefel

Uma comparação da convergência do método de descida do gradiente com tamanho de passo ótimo (em verde) e o método do gradiente conjugado (em vermelho) para a minimização da forma quadrática com um sistema linear dado. O gradiente conjugado, assumindo aritmética exata, converge em no máximo n passos onde n é o tamanho da matriz do sistema (no exemplo, n=2).

Primeiro passo: Escolha $x_0\in {ℝ}^n$
  Se $\mathrm{\nabla }f(x_0)=0$, então pare: $\overline{x}=x_0$
  Senão: $d_0=-\mathrm{\nabla }f(x_0)=-A{x}_0-a$
  Calcular $t_0=\frac{\Vert \mathrm{\nabla }f(x_0){\Vert }^2}{d_0^{\mathrm{\top }}A{d}_0}$
  $x_1=x_0+t_0{d}_0$


Passo iterativo $k$: Dado $x_k\in {ℝ}^n$
  Se $\mathrm{\nabla }f(x_k)=0$, então pare: $\overline{x}=x_k$
  Senão: $d_k=-\mathrm{\nabla }f(x_k)+\frac{\mathrm{\nabla }f(x_k{)}^{\mathrm{\top }}A{d}_k}{d_k^{\mathrm{\top }}A{d}_k}{d}_k$
  $t_k=\frac{\Vert {\mathrm{\nabla }}^{2}f(x_k){\Vert }^2}{d_k^{\mathrm{\top }}A{d}_k}$
  $x_{k+1}=x_k+t_k{d}_k$

Pode-se verificar facilmente que $d_{k+1}\perp d_k$. De fato, como $d_{k+1}=-\mathrm{\nabla }f(x_{k+1})+\frac{\mathrm{\nabla }f(x_{k+1}{)}^{\mathrm{\top }}A{d}_k}{d_k^{\mathrm{\top }}A{d}_k}{d}_k$, tem-se $A{d}_{k+1}=-A\mathrm{\nabla }f(x_{k+1})+\frac{\mathrm{\nabla }f(x_{k+1}{)}^{\mathrm{\top }}A{d}_k}{d_k^{\mathrm{\top }}A{d}_k}A{d}_k$. Logo, $d_k^{\mathrm{\top }}A{d}_{k+1}=-\mathrm{\nabla }f(x_{k+1}{)}^{\mathrm{\top }}A{d}_k+\frac{\mathrm{\nabla }f(x_{k+1}{)}^{\mathrm{\top }}A{d}_k}{d_k^{\mathrm{\top }}A{d}_k}{d}_k^{\mathrm{\top }}A{d}_k=-\mathrm{\nabla }f(x_{k+1}{)}^{\mathrm{\top }}A{d}_k+\mathrm{\nabla }f(x_{k+1}{)}^{\mathrm{\top }}A{d}_k=0$.

---

Exemplos

Considere $f:{ℝ}^2\mapsto ℝ$ definida por $f(x,y)=\frac{1}{2}(x^2+y^2)=\frac{1}{2}\left[\begin{array}{cc}x& y\end{array}\right]\left[\begin{array}{cc}1& 0\\ 0& 1\end{array}\right]\left[\begin{array}{c}x\\ y\end{array}\right]$. Em outros termos, tomando $u=\left[\begin{array}{c}x\\ y\end{array}\right]$, tem-se $f(u)=\frac{1}{2}{u}^{\mathrm{\top }}Au$, onde $A=\left[\begin{array}{cc}1& 0\\ 0& 1\end{array}\right]=I_{2\times 2}$.

Pode-se aplicar o método de direções conjugadas ao seguinte problema

$(P)\{\begin{array}{c}minf(u)\\ u\in {ℝ}^2\end{array}$

Note, desde já, que o conjunto solução é $S=\{\left[\begin{array}{c}0\\ 0\end{array}\right]\}$.

Inicio

• Toma-se $x_0$ arbitrário, por exemplo, $x_0=\left[\begin{array}{c}2\\ 1\end{array}\right]$.
• Avalia-se o gradiente da função $f$ neste ponto inicial: $\mathrm{\nabla }f(x_0)=A{x}_0=I_{2\times 2}{x}_0=x_0$

Iteração 1

• $d_0=-\mathrm{\nabla }f(x_0)=\left[\begin{array}{c}-2\\ -1\end{array}\right]$
• $t_0=\frac{\Vert \mathrm{\nabla }f(x_0){\Vert }^2}{d_0^{\mathrm{\top }}A{d}_0}=\frac{5}{5}=1$
• $x_1=x_0+t_0{d}_0=\left[\begin{array}{c}2\\ 1\end{array}\right]+1\left[\begin{array}{c}-2\\ -1\end{array}\right]=\left[\begin{array}{c}0\\ 0\end{array}\right]$

A seguir, verifica-se se o gradiente se anula no novo ponto $x_1$:

• $\mathrm{\nabla }f(x_1)=A\left[\begin{array}{c}0\\ 0\end{array}\right]=\left[\begin{array}{c}0\\ 0\end{array}\right]$

Como o gradiente já é nulo, não é preciso fazer a segunda iteração, e o ponto $x_1$ é o (único) minimizador global de $f$.

---

Em um caso mais geral, considerando $f:{ℝ}^2\mapsto ℝ$ definida por $f(x,y)=\frac{\lambda }{2}(x^2+y^2)=\frac{1}{2}\left[\begin{array}{cc}x& y\end{array}\right]\left[\begin{array}{cc}\lambda & 0\\ 0& \lambda \end{array}\right]\left[\begin{array}{c}x\\ y\end{array}\right]$, tem-se cálculos muito parecidos em cada passo.

O conjunto solução continua sendo $S=\{\left[\begin{array}{c}0\\ 0\end{array}\right]\}$.

Inicio

• Considere $x_0$ como no primeiro exemplo, ou seja, $x_0=\left[\begin{array}{c}2\\ 1\end{array}\right]$.
• Avalia-se o gradiente da função $f$ neste ponto inicial: $\mathrm{\nabla }f(x_0)=A{x}_0=\lambda x_0$

Iteração 1

• $d_0=-\mathrm{\nabla }f(x_0)=\lambda \left[\begin{array}{c}-2\\ -1\end{array}\right]$
• $t_0=\frac{\Vert \mathrm{\nabla }f(x_0){\Vert }^2}{d_0^{\mathrm{\top }}A{d}_0}=\frac{5{\lambda }^2}{5{\lambda }^3}=\frac{1}{\lambda }$
• $x_1=x_0+\frac{1}{\lambda }\lambda d_0=\left[\begin{array}{c}2\lambda \\ \lambda \end{array}\right]+1\left[\begin{array}{c}-2\lambda \\ -\lambda \end{array}\right]=\left[\begin{array}{c}0\\ 0\end{array}\right]$

A seguir, verifica-se se o gradiente se anula no novo ponto $x_1$:

• $\mathrm{\nabla }f(x_1)=A\left[\begin{array}{c}0\\ 0\end{array}\right]=\lambda \left[\begin{array}{c}0\\ 0\end{array}\right]=\left[\begin{array}{c}0\\ 0\end{array}\right]$

Novamente, o gradiente se anula já na primeira iteração, de modo que $x_1$ é o minimizador global de $f$.

---

Um terceiro exemplo pode ser dado tomando $A=\left[\begin{array}{cc}2& -1\\ -1& 2\end{array}\right]$ e $f:{ℝ}^2\mapsto ℝ$ definida por $f(u)=\frac{1}{2}{u}^{\mathrm{\top }}Au$. Observe que tal matriz é simétrica e definida positiva:

$det(A-\lambda I)=(2-\lambda )(3-\lambda )-1={\lambda }^2-4\lambda -3=(\lambda -3)(\lambda -1)$

Logo, os autovalores de $A$ são $\lambda =1$ e $\lambda =3$. Isso também implica que a função é fortemente convexa.

Aplicando o método:

Início

• Toma-se um ponto arbitrário no plano, por exemplo $x_0=\left[\begin{array}{c}10\\ 20\end{array}\right]$;
• Verifica-se se tal ponto é o minimizador global, avaliando nele o gradiente da função:

$\mathrm{\nabla }f(x_0)=\left[\begin{array}{cc}2& -1\\ -1& 2\end{array}\right]\left[\begin{array}{c}10\\ 20\end{array}\right]=\left[\begin{array}{c}0\\ 30\end{array}\right]≠\left[\begin{array}{c}0\\ 0\end{array}\right]$.

• Já que o gradiente não se anulou no chute inicial, é preciso escolher uma direção e um comprimento de passo para determinar a próxima aproximação:

Iteração 1

$d_0=-\mathrm{\nabla }f(x_0)=\left[\begin{array}{c}0\\ -30\end{array}\right]$

$t_0=\frac{\Vert \left[\begin{array}{cc}0& 30\end{array}\right]{\Vert }^2}{\left[\begin{array}{cc}0& -30\end{array}\right]\left[\begin{array}{cc}2& -1\\ -1& 2\end{array}\right]\left[\begin{array}{c}0\\ -30\end{array}\right]}=\frac{900}{\left[\begin{array}{cc}0& -30\end{array}\right]\left[\begin{array}{c}-30\\ -60\end{array}\right]}=\frac{900}{1800}=\frac{1}{2}$

Feitos esses cálculos, o próximo ponto é dado por

$x_1=x_0+t_0{d}_0=\left[\begin{array}{c}10\\ 20\end{array}\right]+\frac{1}{2}\left[\begin{array}{c}0\\ -30\end{array}\right]=\left[\begin{array}{c}10\\ 20\end{array}\right]+\left[\begin{array}{c}0\\ -15\end{array}\right]=\left[\begin{array}{c}10\\ 5\end{array}\right]$

Para saber se será necessária uma nova iteração, ou se o minimizador foi encontrado, calcula-se o gradiente da função no ponto:

$\mathrm{\nabla }f(x_1)=\left[\begin{array}{cc}2& -1\\ -1& 2\end{array}\right]\left[\begin{array}{c}10\\ 5\end{array}\right]=\left[\begin{array}{c}15\\ 0\end{array}\right]≠\left[\begin{array}{c}0\\ 0\end{array}\right]$.

Novamente, será preciso calcular uma nova direção e um novo comprimento de passo:

Iteração 2

$d_0=\left[\begin{array}{c}-15\\ 0\end{array}\right]+\beta \left[\begin{array}{c}0\\ -30\end{array}\right]=\left[\begin{array}{c}-15\\ -30\beta \end{array}\right]$

onde $\beta$, no algoritmo de Hestenes é dado por:

$\beta =\frac{\left[\begin{array}{cc}15& 0\end{array}\right]\left[\begin{array}{cc}2& -1\\ -1& 2\end{array}\right]\left[\begin{array}{c}0\\ -30\end{array}\right]}{\left[\begin{array}{cc}0& -30\end{array}\right]\left[\begin{array}{cc}2& -1\\ -1& 2\end{array}\right]\left[\begin{array}{c}0\\ -30\end{array}\right]}=\frac{\left[\begin{array}{cc}15& 0\end{array}\right]\left[\begin{array}{c}30\\ -60\end{array}\right]}{\left[\begin{array}{cc}0& -30\end{array}\right]\left[\begin{array}{c}30\\ -60\end{array}\right]}=\frac{15\times 30}{(-30)\times (-60)}=\frac{1}{4}$

Portanto

$d_0=-15\left[\begin{array}{c}1\\ 1/2\end{array}\right]$

Além disso, o tamanho do passo é dado por

$t_1=\frac{\Vert \mathrm{\nabla }f(x_1){\Vert }^2}{d_0^{\mathrm{\top }}A{d}_0}=\frac{15^2}{15^2\left[\begin{array}{cc}1& 1/2\end{array}\right]\left[\begin{array}{cc}2& -1\\ -1& 2\end{array}\right]\left[\begin{array}{c}1\\ 1/2\end{array}\right]}=\frac{1}{\left[\begin{array}{cc}1& 1/2\end{array}\right]\left[\begin{array}{c}3/2\\ 0\end{array}\right]}=\frac{1}{3/2}=\frac{2}{3}$

Portanto

$x_2=x_1+t_1{d}_1=\left[\begin{array}{c}10\\ 5\end{array}\right]-15\frac{2}{3}\left[\begin{array}{c}1\\ 1/2\end{array}\right]=\left[\begin{array}{c}10\\ 5\end{array}\right]-10\left[\begin{array}{c}1\\ 1/2\end{array}\right]=\left[\begin{array}{c}0\\ 0\end{array}\right]$

Obviamente, este é o minimizador procurado (pois o método tem a propriedade de convergência quadrática, ou seja utiliza no máximo $n$ iterações para chegar a solução quando aplicado a funções quadráticas definidas em ${ℝ}^n$)

Implementação em Scilab

Abaixo é apresentada uma implementação deste algoritmo na linguagem de programação utilizada pelo Scilab.

A = [2 1; 1 2];

function [x] = min_gradiente_conjugado(xk)
  //Entrada: xk em R^n, qualquer "chute inicial"
  //  Saída: x, o ponto em que f assume o valor mínimo
  
  k        = 0        //Indica quantos loops já foram feitos
  epsilon  = 0.01
  n        = size(xk,1)
  g        = df(xk)
  seq      = zeros(n,n+1)
  seq(:,1) = xk
  while (norm(g) > epsilon) & (k <= n)
    if (k == 0)
      d = -g
    else
      d = Hestenes(g,d,A)
    end
    t  = busca_linear(g,d,A)
    xk = xk + t*d
    k  = k+1
    seq(:,k+1) = xk
    g  = df(xk)
  end
  plot(seq(1,:),seq(2,:))
  x = xk  
endfunction


function [fu] = f(u)
  fu=(1/2)*(u'*A*u)
endfunction


function [grf] = df(u)
  grf = A*u
endfunction


function [d] = Hestenes(g,d,A)
  d=-g + ((g'*A*d)/(d'*A*d))*d
endfunction


function [t] = busca_linear(g,d,A)
  t=(g'*g)/(d'*A*d)
endfunction

Para facilitar a compreensão do método, pode ser útil exibir as curvas de nível da função. Uma forma de implementar uma função com esse propósito é a seguinte:

function plotar_curvas
  qtd=101
  tam=max(seq)
  x=linspace(-tam,tam,qtd)
  y=x
  z=zeros(qtd,qtd)
  for i=1:qtd
    for j=1:qtd
      z(i,j)=f([x(i);y(j)])
    end
  end
  contour2d(x,y,z,10)
  a=gca();
  a.x_location = "middle"; 
  a.y_location = "middle"; 
endfunction

Algoritmo de Fletcher-Reeves

Esta versão é na verdade uma extensão do algoritmo anterior, permitindo a aplicação no caso de funções que não são quadráticas.

Primeiro passo: Escolha $x_0\in {ℝ}^n$
 Se $\mathrm{\nabla }f(x_0)=0$, então pare: $\overline{x}=x_0$
 Senão: $d_0=-\mathrm{\nabla }f(x_0)$ (como em todo método de descida)
 Calcular $t_0$, através de uma busca linear
 $x_1=x_0+t_0{d}_0$
Passo iterativo:
 Se $\mathrm{\nabla }f(x_k)=0$, então pare: $\overline{x}=x_k$
 Senão: $d_k=-\mathrm{\nabla }f(x_k)+\frac{\Vert \mathrm{\nabla }f(x_k){\Vert }^2}{\Vert \mathrm{\nabla }f(x_{k-1}){\Vert }^2}{d}_{k-1}$
 Calcular $t_k$, através de uma busca linear
 $x_{k+1}=x_k+t_k{d}_k$
 $k=k+1$

Algoritmo de Polak-Ribière

Uma outra versão é a seguinte:

Primeiro passo: Tomar $x_0\in {ℝ}^n$
 Se $\mathrm{\nabla }f(x_0)=0$, então pare: $\overline{x}=x_0$
 Senão: $d_0=-\mathrm{\nabla }f(x_0)$ (como em todo método de descida)
 Calcular $t_0$, através de uma busca linear
 $x_1=x_0+t_0{d}_0$
 $k=1$
Passo iterativo:
 Se $\mathrm{\nabla }f(x_k)=0$, então pare: $\overline{x}=x_k$
 Senão: $d_k=-\mathrm{\nabla }f(x_k)+\frac{\mathrm{\nabla }f(x_k{)}^{\mathrm{\top }}(\mathrm{\nabla }f(x_k)-\mathrm{\nabla }f(x_{k-1}))}{\Vert \mathrm{\nabla }f(x_{k-1}){\Vert }^2}{d}_{k-1}$
 Calcular $t_k$, através de uma busca linear
 $x_{k+1}=x_k+t_k{d}_k$
 $k=k+1$

Algoritmo auxiliar

Para o caso de funções não quadráticas, é preciso usar algum método de busca linear para a implementação do método dos gradientes conjugados, seja a versão de Fletcher-Reeves ou a de Polak-Ribière. Uma possibilidade é a busca de linear de Armijo (ver Izmailov & Solodov (2007), vol 2, pag. 65), cujo algoritmo é esboçado a seguir:

function busca_linear_Armijo (f, theta, alpha, delta, t0)
  while (alpha * pred > ared)
    t = d * t
  end
endfunction

com:

• $pred=-t\theta$
• $\theta (t)=f(x+td)$
• ${\theta }^{\prime }(t)=\mathrm{\nabla }f(x+td{)}^{\mathrm{\top }}d$