Sensitivity of LM minimiser to 'physical' choices

[Ian Culverwell]

The 'physical' limitations on the state vector in the differenced bending angle LM solver, namely

    DO i = 1, nstate

      x%state(i) = x%state(i) + SIGN( MIN( ABS(delta_x(i)), state_sigma(i) ), delta_x(i) ) 

      SELECT CASE( 1 + ((i-1)/nlayer) )

        CASE( 1 )
          IF ( x%state(i) < ropp_ZERO ) THEN
            CALL message(msg_warn, "Levenberg-Marquardt solver returns " // &
              "ne_peak < 0 ... resetting to 1% of background error. \n")
            x%state(i) = 0.01_wp * state_sigma(i)
          END IF

        CASE( 2 )
          IF ( x%state(i) < 0.01_wp*state_sigma(i) ) THEN
            CALL message(msg_info, "Levenberg-Marquardt solver returns " // &
              "r_peak < 10% of background error ... resetting to 10% of background error. \n")
            x%state(i) = 0.10_wp * state_sigma(i)
          END IF

        CASE( 3 )
          IF ( x%state(i) < 0.10_wp*state_sigma(i) ) THEN
            CALL message(msg_info, "Levenberg-Marquardt solver returns " // &
              "h_zero < 10% of background error ... resetting to 10% of background error. \n")
            x%state(i) = 0.10_wp * state_sigma(i)
          END IF

        CASE( 4 )
          IF ( x%state(i) < 1.0E-10_wp*state_sigma(i) ) THEN
            CALL message(msg_info, "Levenberg-Marquardt solver returns " // &
              "h_grad < 1e-10*background error ... resetting to 1e-10*background error. \n")
            x%state(i) = 1.0E-10_wp * state_sigma(i)
          END IF

      END SELECT

    ENDDO

are open to question, to say the least. For example, Sean has had some success with

    x%state(i) = x%state(i) + SIGN( MIN( ABS(delta_x(i)), 0.75_wp*ABS(x%state(i)) ), delta_x(i) )

These choices need careful consideration and experimentation. This ticket is the place to hold the results.

[Ian Culverwell]

Note the typo in the code snippet above. Instead of

        CASE( 2 )
          IF ( x%state(i) < 0.01_wp*state_sigma(i) ) THEN
            CALL message(msg_info, "Levenberg-Marquardt solver returns " // &
              "r_peak < 10% of background error ... resetting to 10% of background error. \n")
            x%state(i) = 0.10_wp * state_sigma(i)
          END IF

it should read

         CASE( 2 )
         IF ( x%state(i) < 0.10_wp*state_sigma(i) ) THEN
            CALL message(msg_info, "Levenberg-Marquardt solver returns " // &
              "r_peak < 10% of background error ... resetting to 10% of background error. \n")
            x%state(i) = 0.10_wp * state_sigma(i)
          END IF

(i.e. 10% of sigma in the IF statement, not 1%.) This bug has been fixed in the ROPP-11 prototype at r6903.

However, x%state(2) (etc) is the radial distance to the peak of the VaryChap layer, say around 6400 + 300 = 6700 km, while state_sigma(2) (etc) is around 50-150 km, and therefore neither condition is ever met in practice. It might be better to replace this condition by

          IF ( x%state(i) - offset(i) < 0.10_wp*state_sigma(i) ) THEN
            CALL message(msg_info, "Levenbetimerg-Marquardt solver returns " // &
              "r_peak - offset < 10% of background error ... resetting to 10% of background error. \n")
            x%state(i) = offset(i) + 0.10_wp * state_sigma(i)
          END IF

where offset(i) is something like R_earth or RoC or the r_peak of the background. The value of offset would depend on whether preconditioning is on or off. (A quick test of this, using offset=6.4E6_wp, led to an improvement in the 1 layer retrieval of the Metop-A extension dataset GRAS_1B_M02_20200801185622Z_20200801190137Z_N_T_20210507160639Z_G11_NN_iono.nc, but that was the only occasion (out of 2x19 = 38 retrievals) that the check was invoked. So it's unlikely to revolutionise the results.)

[Ian Culverwell]

For the record, these are the results of using the old code:

And these are the results of using the new code:

[Ian Culverwell]

G11_1L_old.png

[Ian Culverwell]

G11_1L_new.png