access to last hessian in BoundedBroydenFletcherGoldfarbShanno

[aeternam7]

hi, i'm using BoundedBroydenFletcherGoldfarbShanno optimization class for modelling GARCH(1,1)

i get the parameter estimates of the model, but also need the standard errors of these estimates

these standard errors are computed with the hessian matrix in the convergence point

i do not add any hessian definition, so i suppose it is computed numerically

so, how can i access to last hessian computed in BoundedBroydenFletcherGoldfarbShanno class?

thanks in advance

[cesarsouza]

Hi @aeternam7,

The BoundedBroydenFletcherGoldfarbShannon class implements a limited-memory bounded BFGS (L-BFGS-b) algorithm, meaning the Hessian is never completely computed and stored in memory during its execution. As such, it is not really possible to access the Hessian matrix from BBFGS at this time.

It may be possible to compute the Hessian matrix outside of the algorithm though. I am trying to determine how this could be done, but if you find something, please update the issue as well!

Regards,
Cesar

[aeternam7]

hi, cesar, thanks for your reply

i think i'll translate this R code to vb.net

https://rdrr.io/cran/CDM/src/R/numerical_Hessian.R

it gives the same result as the referece applications that i'm using to test my code

regards

[cesarsouza]

Hi @aeternam7,

I was about to add a new Hessian(x) method to the FiniteDifferences class. Would you like to test it?

I will commit it in a few minutes.

Regards,
Cesar

[aeternam7]

of course, i'll test it and compare with the results of some statistical packages

regards

[cesarsouza]

I've just added, but forgot to reference the issue in the commit message. The new method can be seen here:

framework/Sources/Accord.Math/Differentiation/FiniteDifferences.cs

Line 287 in 3131c41

public double[][] Hessian(double[] x)

If you don't want to recompile the entire framework, feel free to just copy & paste that file in your project with a different name such as 'FineDifferences2' (so it does not conflict with the current FiniteDifferences in Accord.Math).

Regards,
Cesar

[aeternam7]

thanks

i tested your code

the result for the case i'm applying is

hess [0][0]= 206071,54949922
hess [0][1]= 788196,943311202
hess [0][2]= -79,6849089667598
hess [0][3]= 2499,39552610101
hess [1][0]= 788196,943311202
hess [1][1]= 3902120617,83357
hess [1][2]= 228359,154267735
hess [1][3]= 836441,535323473
hess [2][0]= -79,6849089667598
hess [2][1]= 228359,154267735
hess [2][2]= 569,324365336775
hess [2][3]= 606,761045782167
hess [3][0]= 2499,39552125697
hess [3][1]= 836441,535316759
hess [3][2]= 606,761045771678
hess [3][3]= 86609609,5954028

its quite different than the result of this code

https://rdrr.io/cran/CDM/src/R/numerical_Hessian.R

        [,1]         [,2]        [,3]        [,4]
        [1,]  289452.471    7575924.6   -419.4690   2716.6562
        [2,] 7575924.618 2895928593.8 210342.3287 764151.2843
        [3,]    -419.469     210342.3    569.0492    605.1528
        [4,]    2716.656     764151.3    605.1528    977.6352

i translated the above code to vb; so you can see how it works

   `Private Function fNumericalHessianAccord(par As Double(), fun As Func(Of Double(), Double), Optional h As Double = 0.00001, Optional bGradient As Boolean = True) As Double()()

    Dim NP As Integer = par.Length
    Dim f0 As Double = fun(par)

    Dim fh(NP - 1) As Double
    For i = 0 To NP - 1
        fh(i) = Double.NaN
    Next

    Dim f2h(NP - 1) As Double
    fh.CopyTo(f2h, 0)

    Dim hess(NP - 1)() As Double
    For i = 0 To NP - 1
        Dim hessRow(NP - 1) As Double
        fh.CopyTo(hessRow, 0)
        hess(i) = hessRow
    Next

    Dim fhh(NP - 1)() As Double
    For i = 0 To NP - 1
        Dim fhhRow(NP - 1) As Double
        fh.CopyTo(fhhRow, 0)
        fhh(i) = fhhRow
    Next

    'select h parameters according to size of parameters
    Dim abs_par(NP - 1) As Double
    For i = 0 To NP - 1
        abs_par(i) = Math.Abs(par(i))
    Next

    Dim hvec(NP - 1) As Double
    For i = 0 To NP - 1
        If abs_par(i) > 1 Then
            hvec(i) = h * abs_par(i)
        Else
            hvec(i) = h
        End If
    Next


    'loop for computing f(x+h)
    Dim par1(NP - 1) As Double
    For ii = 0 To NP - 1
        par.CopyTo(par1, 0)
        par1(ii) = par(ii) + hvec(ii)
        fh(ii) = fun(par1)
    Next

    'computation of the gradient
    Dim grad1(NP - 1) As Double
    If bGradient Then
        For ii = 0 To NP - 1
            grad1(ii) = (fh(ii) - f0) / hvec(ii)
        Next
    End If

    'second partial derivatives
    Dim fh1(NP - 1)() As Double
    For i = 0 To NP - 1
        Dim fh1Row(NP - 1) As Double
        fh.CopyTo(fh1Row, 0)
        fh1(i) = fh1Row
    Next

    'igual pero la transpuesta
    Dim fh2(NP - 1)() As Double
    For i = 0 To NP - 1
        Dim fh2Row(NP - 1) As Double
        For j = 0 To NP - 1
            fh2Row(j) = fh(i)
        Next
        fh2(i) = fh2Row
    Next

    ' computation f(x+2*h)
    For ii = 0 To NP - 1
        par.CopyTo(par1, 0)
        par1(ii) = par(ii) + 2 * hvec(ii)
        f2h(ii) = fun(par1)
    Next

    'computation f(x+h,y+h) 
    For ii = 0 To NP - 1
        For jj = 0 To NP - 1
            If (ii < jj) Then
                par.CopyTo(par1, 0)
                par1(ii) = par(ii) + hvec(ii)
                par1(jj) = par(jj) + hvec(jj)
                fhh(jj)(ii) = fun(par1)
                fhh(ii)(jj) = fhh(jj)(ii)
            End If
        Next
    Next

    'hacemos el outer product de los vectores hvec
    Dim h_squared(NP - 1)() As Double
    For ii = 0 To NP - 1
        Dim h_squaredRow(NP - 1) As Double
        For jj = 0 To NP - 1
            If (ii <= jj) Then 'diagonal y bajo la diagonal, y se replica luego sobre la diagonal
                h_squaredRow(jj) = hvec(ii) * hvec(jj)
            End If
        Next
        h_squared(ii) = h_squaredRow
    Next
    'replicamos sobre la diagonal
    For ii = 0 To NP - 1
        For jj = 0 To NP - 1
            If (ii > jj) Then
                h_squared(ii)(jj) = h_squared(jj)(ii)
            End If
        Next
    Next


    'rellenamos solo la diagonal superior, y replicamos abajo
    For ii = 0 To NP - 1
        For jj = 0 To NP - 1
            If (ii < jj) Then
                hess(ii)(jj) = (fhh(ii)(jj) - fh1(ii)(jj) - fh2(ii)(jj) + f0) / h_squared(ii)(jj)
                hess(jj)(ii) = hess(ii)(jj)
            End If
        Next
    Next
    'ahora hacemos la diagonal del hessiano
    For ii = 0 To NP - 1
        hess(ii)(ii) = (f2h(ii) - 2 * fh(ii) + f0) / hvec(ii) ^ 2
    Next

    Return hess

End Function`

fGarch package for R show these results

              mu         omega       alpha1        beta1
 mu     -2.888344e+05   -11132706.2     419.9592   -2730.8111
 omega  -1.113271e+07 -7312863526.1 -252050.0002 -927298.2752
 alpha1  4.199592e+02     -252050.0    -569.0740    -605.1578
 beta1  -2.730811e+03     -927298.3    -605.1578    -977.7221

but i still don't know which of these functions is using this package to compute the hessian

.garchRoptimhess <-
function(par, .params, .series, eps = 1.0e-4)
{
# A function implemeted by Diethelm Wuertz

# Description:
#   Compute Hessian via R's function optimHess()

# Arguments:
#   par -
#   .params -
#   .series -
#   eps -

# FUNCTION:

# Take Time:
.StartHessian <- Sys.time()

# Compute Hessian:
H <- optimHess(par, .garchLLH)
H <- 0.5 * (H + t(H))
nm <- names(par)
dimnames(H) <- list(nm, nm)

# Elapsed Time:
time = Sys.time() - .StartHessian
attr(H, "time") = time

# Return Value:
H
}

.garchRCDAHessian <-
function(par, .params, .series, eps = 1.0e-4)
{
# A function implemented by Diethelm Wuertz

# Description:
#   Compute CDA (central difference approximated) Hessian

# Arguments:
#   par -
#   .params -
#   .series -
#   eps -

# Reference:
#   http://status.sph.umich.edu/computing/manuals/sas8/stat/chap46/sect26.htm

# FUNCTION:

# Take start time
.StartHessian <- Sys.time()

# Algorithm:
algorithm = .params$control$algorithm[1]
.trace = FALSE

# Compute Hessian:
eps = eps * par
n = length(par)
H = matrix(0, ncol = n, nrow = n)
for (i in 1:n) {
    for (j in 1:n) {
        x1 = x2 = x3 = x4 = par
        x1[i] = x1[i] + eps[i]
        x1[j] = x1[j] + eps[j]
        x2[i] = x2[i] + eps[i]
        x2[j] = x2[j] - eps[j]
        x3[i] = x3[i] - eps[i]
        x3[j] = x3[j] + eps[j]
        x4[i] = x4[i] - eps[i]
        x4[j] = x4[j] - eps[j]
        H[i, j] = (
         .garchLLH(x1, .trace) -
         .garchLLH(x2, .trace) -
         .garchLLH(x3, .trace) +
         .garchLLH(x4, .trace) ) / (4*eps[i]*eps[j])
    }
}
colnames(H) = rownames(H) = names(par)

# Attribute execution time:
time = Sys.time() - .StartHessian
attr(H, "time") = time

# Return Value:
H
 }

 .garchTSHessian <-
function(par, .params, .series, eps = NA)
{
# A function implemented by Diethelm Wuertz

# Description:
#   Compute two sided (TS) approximated Hessian

# Arguments:
#   par -
#   .params -
#   .series -
#   eps - not used

# FUNCTION:

# Take start time
.StartHessian <- Sys.time()

# Algorithm:
algorithm = .params$control$algorithm[1]

# Compute Hessian:
H <- .hessian2sided(f = .garchLLH, x = par, trace = FALSE, fGarchEnv = FALSE)
colnames(H) = rownames(H) = names(par)

# Attribute execution time:
time = Sys.time() - .StartHessian
attr(H, "time") = time

# Return Value:
H
}

.hessian2sided <-
function(f, x, ...)
{
# A function adapted from Kevin Sheppard's Matlab garch toolbox
# ... implemented by Alexios Ghalanos in his rgarch package
# ... R port for Rmetrics' fGarch by Diethelm Wuertz

# Description:
#   Computes two sided (TS) approximated Hessian

# Arguments:
#   f -
#   x -

# Notes:
#   requires package Matrix (added as suggestion)

# FUNCTION:

# Settings:
n = length(x)
fx <- f(x, ...)
eps = .Machine$double.eps

# Compute the stepsize (h)
h = eps^(1/3) *
    apply( as.data.frame(x), 1, FUN = function(z) max(abs(z), 1.0e-2))
xh = x + h
h = xh - x
ee = Matrix::Matrix(diag(h), sparse = TRUE)

# Compute forward and backward steps:
gp = vector(mode = "numeric", length = n)
for(i in 1:n) gp[i] <- f(x + ee[, i], ...)
gm = vector(mode = "numeric", length = n)
for(i in 1:n) gm[i] <- f(x - ee[, i], ...)
H = h %*% t(h)
Hm = H
Hp = H

# Compute "double" forward and backward steps:
for(i in 1:n){
    for(j in  i:n){
        Hp[i, j] <- f(x + ee[, i] + ee[, j], ...)
        Hp[j, i] <- Hp[i, j]
        Hm[i, j] <- f(x - ee[, i] - ee[, j], ...)
        Hm[j, i] <- Hm[i, j]
    }
}

# Compute the hessian:
for(i in 1:n){
    for(j in  i:n){
        H[i, j] = ((Hp[i, j] - gp[i] - gp[j] + fx + fx - gm[i] - gm[j] +
                    Hm[i, j]) / H[i, j]) / 2
        H[j, i] = H[i, j]
    }
}

# Return Value:
return(H)
}

i cant test your code against other packages, but...
sincerely, i have not the knowledgment to determine which algorithm is better
if you tell me your code is better, i will use it :)

[cesarsouza]

Thanks a lot @aeternam7!

I am also still testing it, but can you try to increase the stepSize parameter? By default the code is using 1e-2, which might be quite low compared to the 0.00001 in the VB version.

i.e. you can use

var calculator = new FiniteDifferences(variables: 4, function: function, order: 1, stepSize: 1e-5);

Also, do you have an example I can use to reproduce the numbers you have shown above?

Regards,
Cesar

[cesarsouza]

For example, do you think you could provide the arguments you have passed to the parameters of the functions from R's fGarch package?

[aeternam7]

Hi, Cesar, fGarch is not suitable for the test, because i don’t know exactly which formula is trying to maximize; it’s a recursive formula for maximum likelihood that depends on initial values; i also don’t know which initial values is using fGarch; i get similar results as fGarch but don’t know if the differences are because of the formula, because of the initial values, or because of the optimization algorithm results; so it should not be used as a test reference

If you are interested in the fGarh formula, it is on fGarch/src/llhGarch.f. It depends on many parameters, as fGarch computes a wide range of different GARCH models

Anyway, tell me which of the R hessian methods i showed above you are interested in for the test; and i will translate it to VB.NET

[cesarsouza]

Hi @aeternam7,

What I mean is, can you please provide an example of the arguments that you have passed to at least one of the examples you have posted above? It can be, for example, the arguments you have passed to R's https://rdrr.io/cran/CDM/src/R/numerical_Hessian.R. I just want to be able to call a code example that can give me results that are similar to the ones you would have expected, like for example the ones you have shown in

           [,1]         [,2]        [,3]        [,4]
[1,]  289452.471    7575924.6   -419.4690   2716.6562
[2,] 7575924.618 2895928593.8 210342.3287 764151.2843
[3,]    -419.469     210342.3    569.0492    605.1528
[4,]    2716.656     764151.3    605.1528    977.6352

Regards,
Cesar

[aeternam7]

Hi, i passed a recursive formula which operates with an array of 1.300 observations … so lets create a new better test

code with R hessian

 `     Module Module1

  Sub main()

    Dim aParam As Double() = {1.790978, 0.00009408872, 0.9888748, 0.00000001}
    Dim hess()() As Double = fNumericalHessianAccord(aParam, funcionObjetivoLBFGS)

End Sub


Private Function fNumericalHessianAccord(par As Double(), fun As Func(Of Double(), Double), Optional h As Double = 0.00001, Optional bGradient As Boolean = True) As Double()()

    Dim NP As Integer = par.Length
    Dim f0 As Double = fun(par)

    Dim fh(NP - 1) As Double
    For i = 0 To NP - 1
        fh(i) = Double.NaN
    Next

    Dim f2h(NP - 1) As Double
    fh.CopyTo(f2h, 0)

    Dim hess(NP - 1)() As Double
    For i = 0 To NP - 1
        Dim hessRow(NP - 1) As Double
        fh.CopyTo(hessRow, 0)
        hess(i) = hessRow
    Next

    Dim fhh(NP - 1)() As Double
    For i = 0 To NP - 1
        Dim fhhRow(NP - 1) As Double
        fh.CopyTo(fhhRow, 0)
        fhh(i) = fhhRow
    Next

    'select h parameters according to size of parameters
    Dim abs_par(NP - 1) As Double
    For i = 0 To NP - 1
        abs_par(i) = Math.Abs(par(i))
    Next

    Dim hvec(NP - 1) As Double
    For i = 0 To NP - 1
        If abs_par(i) > 1 Then
            hvec(i) = h * abs_par(i)
        Else
            hvec(i) = h
        End If
    Next


    'loop for computing f(x+h)
    Dim par1(NP - 1) As Double
    For ii = 0 To NP - 1
        par.CopyTo(par1, 0)
        par1(ii) = par(ii) + hvec(ii)
        fh(ii) = fun(par1)
    Next

    'computation of the gradient
    Dim grad1(NP - 1) As Double
    If bGradient Then
        For ii = 0 To NP - 1
            grad1(ii) = (fh(ii) - f0) / hvec(ii)
        Next
    End If

    'second partial derivatives
    Dim fh1(NP - 1)() As Double
    For i = 0 To NP - 1
        Dim fh1Row(NP - 1) As Double
        fh.CopyTo(fh1Row, 0)
        fh1(i) = fh1Row
    Next

    'igual pero la transpuesta
    Dim fh2(NP - 1)() As Double
    For i = 0 To NP - 1
        Dim fh2Row(NP - 1) As Double
        For j = 0 To NP - 1
            fh2Row(j) = fh(i)
        Next
        fh2(i) = fh2Row
    Next

    ' computation f(x+2*h)
    For ii = 0 To NP - 1
        par.CopyTo(par1, 0)
        par1(ii) = par(ii) + 2 * hvec(ii)
        f2h(ii) = fun(par1)
    Next

    'computation f(x+h,y+h) 
    For ii = 0 To NP - 1
        For jj = 0 To NP - 1
            If (ii < jj) Then
                par.CopyTo(par1, 0)
                par1(ii) = par(ii) + hvec(ii)
                par1(jj) = par(jj) + hvec(jj)
                fhh(jj)(ii) = fun(par1)
                fhh(ii)(jj) = fhh(jj)(ii)
            End If
        Next
    Next

    'hacemos el outer product de los vectores hvec
    Dim h_squared(NP - 1)() As Double
    For ii = 0 To NP - 1
        Dim h_squaredRow(NP - 1) As Double
        For jj = 0 To NP - 1
            If (ii <= jj) Then 'diagona y bajo la diagonal, y se replica para sobre diagonal
                h_squaredRow(jj) = hvec(ii) * hvec(jj)
            End If
        Next
        h_squared(ii) = h_squaredRow
    Next
    'replicamos a la diagonal superior
    For ii = 0 To NP - 1
        For jj = 0 To NP - 1
            If (ii > jj) Then
                h_squared(ii)(jj) = h_squared(jj)(ii)
            End If
        Next
    Next


    'rellenamos solo la diagonal superior, y replicamos abajo
    For ii = 0 To NP - 1
        For jj = 0 To NP - 1
            If (ii < jj) Then
                hess(ii)(jj) = (fhh(ii)(jj) - fh1(ii)(jj) - fh2(ii)(jj) + f0) / h_squared(ii)(jj)
                hess(jj)(ii) = hess(ii)(jj)
            End If
        Next
    Next
    'ahora hacemos la diagonal del hessiano
    For ii = 0 To NP - 1
        hess(ii)(ii) = (f2h(ii) - 2 * fh(ii) + f0) / hvec(ii) ^ 2
    Next


    For i = 0 To 3
        For j = 0 To 3
            Debug.WriteLine("hess [" & i & "][" & j & "]= " & hess(i)(j))
        Next
    Next



    Return hess

End Function



Dim funcionObjetivoLBFGS As Func(Of Double(), Double) = Function(x)

                                                            Dim iSizeArray = 10
                                                            Dim arrayDecimalesTrabajo() As Double = {2.2227, 2.2188, 2.2144, 2.204, 2.2006, 2.2053, 2.2053, 2.2035, 2.1969, 2.2033}

                                                            Dim mu As Double = x(0)
                                                            Dim omega As Double = x(1)
                                                            Dim alpha As Double = x(2)
                                                            Dim beta As Double = x(3)

                                                            Dim e2(iSizeArray - 1) As Double
                                                            Dim mediae2 As Double = 0.0
                                                            For i = 0 To iSizeArray - 1
                                                                e2(i) = (arrayDecimalesTrabajo(i) - mu) ^ 2
                                                                mediae2 = mediae2 + e2(i)
                                                            Next

                                                            mediae2 = mediae2 / iSizeArray
                                                            Dim sigmat2(iSizeArray - 1) As Double
                                                            sigmat2(0) = omega + alpha * mediae2
                                                            For i = 1 To iSizeArray - 1
                                                                sigmat2(i) = omega + alpha * e2(i - 1)
                                                            Next

                                                            sigmat2(0) = sigmat2(0) + beta * mediae2
                                                            For i = 1 To iSizeArray - 1
                                                                sigmat2(i) = sigmat2(i) + beta * sigmat2(i - 1)
                                                            Next

                                                            Dim sumalikelihood As Double = iSizeArray * Math.Log(2 * Math.PI)
                                                            For i = 0 To iSizeArray - 1
                                                                Dim dMathlog As Double = Math.Log(sigmat2(i))
                                                                If Double.IsInfinity(dMathlog) Then Exit For
                                                                If Double.IsNaN(dMathlog) Then Exit For
                                                                sumalikelihood = sumalikelihood + dMathlog + e2(i) / sigmat2(i)
                                                            Next
                                                            Return 0.5 * sumalikelihood
                                                        End Function

 End Module`

code with accord, c#

` class Program
{

    static void Main(string[] args)
    {
      double[] aParam = {1.790978, 9.408872E-05, 0.9888748, 1E-08};
        

        Func<double[], double> funcionObjetivoLBFGS = x =>
        {
            int iSizeArray = 10;
            double[] arrayDecimalesTrabajo = { 2.2227, 2.2188, 2.2144, 2.204, 2.2006, 2.2053, 2.2053, 2.2035, 2.1969, 2.2033 };
            
            double mu = x[0];
            double omega = x[1];
            double alpha = x[2];
            double beta = x[3];
            
            double[] e2 = new double[iSizeArray];
              double mediae2 = 0.0;
            for (int i = 0; i <= iSizeArray - 1; i++)
            {
                e2[i] = Math.Pow((arrayDecimalesTrabajo[i] - mu), 2);
                mediae2 = mediae2 + e2[i];
            }

            mediae2 = mediae2 / iSizeArray;
            double[] sigmat2 = new double[iSizeArray];
            sigmat2[0] = omega + alpha * mediae2;
            for (int i = 1; i <= iSizeArray - 1; i++)
            {
                sigmat2[i] = omega + alpha * e2[i - 1];
            }
            sigmat2[0] = sigmat2[0] + beta * mediae2;
            for (int i = 1; i <= iSizeArray - 1; i++)
            {
                sigmat2[i] = sigmat2[i] + beta * sigmat2[i - 1];
            }

            double sumalikelihood = iSizeArray * Math.Log(2 * Math.PI);
            for (int i = 0; i <= iSizeArray - 1; i++)
            {
                double dMathlog = Math.Log(sigmat2[i]);
                if (double.IsInfinity(dMathlog))
                    break; 
                if (double.IsNaN(dMathlog))
                    break; 
                sumalikelihood = sumalikelihood + dMathlog + e2[i] / sigmat2[i];
            }
            return 0.5 * sumalikelihood;

        };
        
        var calculator = new FiniteDifferences2(variables: 4, function: funcionObjetivoLBFGS, order: 1, stepSize: 1e-5);
        
        double[] result = calculator.Compute(aParam);
        double[][] hess = calculator.Hessian(aParam);


        for (int i = 0; i <= 3; i++)
        {
            for (int j = 0; j <= 3; j++)
            {
                Debug.WriteLine("hess [" + i + "][" + j + "]= " +  hess[i][j]);
            }
         
        }
        
        System.Console.ReadKey();
        
    }
            
}`

[aeternam7]

sorry, the results

for vb.net

hess [0][0]= -57,6160442364672
hess [0][1]= -1,27840745168738
hess [0][2]= 0,0255249049748045
hess [0][3]= 0,0231494572415752
hess [1][0]= -1,27840745168738
hess [1][1]= 173,045435758468
hess [1][2]= 29,9809421733244
hess [1][3]= 29,8568458845239
hess [2][0]= 0,0255249049748045
hess [2][1]= 29,9809421733244
hess [2][2]= 5,2012616436059
hess [2][3]= 5,17856868498256
hess [3][0]= 0,0231494572415752
hess [3][1]= 29,8568458845239
hess [3][2]= 5,17856868498256
hess [3][3]= 5,15878895157584

for accord

hess [0][0]= -57,6110556635488
hess [0][1]= -1,2913340264809
hess [0][2]= 0,0255249620040496
hess [0][3]= 247,961651211881
hess [1][0]= -1,29133371408033
hess [1][1]= 751,940878763617
hess [1][2]= 29,9744342530767
hess [1][3]= -4719898,51563714
hess [2][0]= 0,0255249621675288
hess [2][1]= 29,9744342523859
hess [2][2]= 5,20144178726767
hess [2][3]= 0
hess [3][0]= 247,959053572999
hess [3][1]= -4719898,51546564
hess [3][2]= 0
hess [3][3]= -44404036008,1383