docs: private void computeDandH( DenseMatrix64F param , DenseMatrix64F x , DenseMatrix64F y )
      {
          func.compute(param,x, tempDH);
          subtractEquals(tempDH, y);
  
          if (jacobianFactory != null)
          {
              jacobianFactory.computeJacobian(param, x, jacobian);
          }
          else
          {
              computeNumericalJacobian(param,x,jacobian);
          }
  
          int numParam = param.getNumElements();
          int length = y.getNumElements();
  
          // d = average{ (f(x_i;p) - y_i) * jacobian(:,i) }
          for( int i = 0; i < numParam; i++ ) {
              double total = 0;
              for( int j = 0; j < length; j++ ) {
                  total += tempDH.get(j,0)*jacobian.get(i,j);
              }
              d.set(i,0,total/length);
          }
  
          // compute the approximation of the hessian
          multTransB(jacobian,jacobian,H);
          scale(1.0/length,H);
      }