code: def line_search_armijo(f, xk, pk, gfk, old_fval,
                         args=(), c1=1e-4, alpha0=0.99):
      
      xk = np.atleast_1d(xk)
      fc = [0]
  
      def phi(alpha1):
          fc[0] += 1
          return f(xk + alpha1 * pk, *args)
  
      if old_fval is None:
          phi0 = phi(0.)
      else:
          phi0 = old_fval
  
      derphi0 = np.sum(pk * gfk)  # Quickfix for matrices
      alpha, phi1 = scalar_search_armijo(
          phi, phi0, derphi0, c1=c1, alpha0=alpha0)
  
      return alpha, fc[0], phi1
docs: Armijo linesearch function that works with matrices
  
      find an approximate minimum of f(xk+alpha*pk) that satifies the
      armijo conditions.
  
      Parameters
      ----------
  
      f : function
          loss function
      xk : np.ndarray
          initial position
      pk : np.ndarray
          descent direction
      gfk : np.ndarray
          gradient of f at xk
      old_fval : float
          loss value at xk
      args : tuple, optional
          arguments given to f
      c1 : float, optional
          c1 const in armijo rule (>0)
      alpha0 : float, optional
          initial step (>0)
  
      Returns
      -------
      alpha : float
          step that satisfy armijo conditions
      fc : int
          nb of function call
      fa : float
          loss value at step alpha