code: def _laplace_approximation_newton(self, tol=1e-6, stepsz=0.9, verbose=False):
          
          from pylds.util import solve_symm_block_tridiag, scipy_solve_symm_block_tridiag
          scale = self.T * self.D_emission
  
          def newton_step(x, stepsz):
              assert 0 <= stepsz <= 1
              g = self.gradient_log_joint(x)
              H_diag, H_upper_diag = self.sparse_hessian_log_joint(x)
              Hinv_g = -scipy_solve_symm_block_tridiag(-H_diag / scale,
                                                 -H_upper_diag / scale,
                                                 g / scale)
              return x - stepsz * Hinv_g
  
          if verbose:
              print("Fitting Laplace approximation")
  
          itr = [0]
          def cbk(x):
              print("Iteration: ", itr[0],
                    "\tObjective: ", (self.log_joint(x) / scale).round(4),
                    "\tAvg Grad: ", (self.gradient_log_joint(x).mean() / scale).round(4))
              itr[0] += 1
  
          # Solve for optimal x with Newton's method
          x = self.gaussian_states
          dx = np.inf
          while dx >= tol:
              xnew = newton_step(x, stepsz)
              dx = np.mean(abs(xnew - x))
              x = xnew
  
              if verbose:
                  cbk(x)
  
          assert np.all(np.isfinite(x))
          if verbose:
              print("Done")
  
          return x
docs: Solve a block tridiagonal system with message passing.