ralna / RALFit

A non-linear least squares solver
Other
25 stars 6 forks source link

Recursive call with regularization (Tensor Newton) #87

Open talassio opened 4 years ago

talassio commented 4 years ago

Tensor Newton subproblem solve (with regularization)

Part 1 of 2 Nan's under valid options in recursive call (regularization)

In the recursive call (model=4) the starting point is either x=d(:)=0 or x=d(:)=1.0e-12, in the zero case, we have that update_regularized_gradient(g,x,normx,options) will have normx=0 and

...
Case (2)
             g = g - options%regularization_term*(normx**(options%                &
            regularization_power-2.0_wp))*x

will produce NaN when 0^0 is encountered (which is a default valid option combination...)

Fix (pushed to box-logic)

     Select Case (options%regularization)
      Case (1)
        g = g - options%regularization_term*x
      Case (2)
        If (normx==0.0_wp) Then
!         Avoid undefined behaviour (Knuth 1992, 0^0 := 1)
          g = g - options%regularization_term*x
        Else
          g = g - options%regularization_term*(normx**(options%                &
            regularization_power-2.0_wp))*x
        End If
      End Select

This is pushed to box-logic.

Part 2 of 2 Subsolvers fail under default starting point

Further to this problem, d(:)=0.0 seems not to be a valid starting point for the sub-solvers

      0 4.9210E+01  4.43299E+01  4.46844E+00  1.00E+02 -1.00E+00 --I-
*** Solving the regularized subproblem using Galahad DRQS ***
*** Subproblem solution not found.  Trying next method ***
*** Solving the regularized subproblem using RALFit Solver ***
*** Error or Subproblem solution NOT found ***
 -------------------------------------------------------------------------------
 Status: aborted, no progress

yet forcing x=d(:) = 1e-12 won't make the subsolvers fail but it neither provides any progress in the solve (need to look into)

Potentially, this can be tackled by using a favourable starting point for the recursive solve e.g. x=d(:)=sqrt(norm(w%g))??? (probably not, need to think about it)

For now the solve just fails with no progress. Options used:

!     Set options
      Call e04zmf(handle,'BXNL Use Second Derivatives = Yes',ifail)
      Call e04zmf(handle,'BXNL Model = Tensor-Newton',ifail)
      Call e04zmf(handle,'BXNL NLLS Method = AINT',ifail)
      Call e04zmf(handle,'BXNL Glob Method = REG',ifail)
      Call e04zmf(handle,'Bxnl Basereg Type= EXPAND-1-DOF',ifail)
      Call e04zmf(handle,'Bxnl TN METHOD= Implicit',ifail)
talassio commented 4 years ago

Part 1 just committed to box-logic branch.