Minimizing the Cauchy function

A good trick (I discovered accidentally) is to use the weight  

$$\bold W \eq \ {\bf diag} \left( {1\over{\sqrt{1+r_i^2/\bar r^2}}} \right)$$ (16)

Sergey Fomel points out that this weight arises from minimizing the Cauchy function:  

$$E(\bold r) \eq \sum_i \ \log (1+r_i^2/\bar r^2)$$ (17)

A plot of this function is found in Figure 3.

 

cauchy
Figure 3 The coordinate is m. Top is Cauchy measures of m-1. Bottom is the same measures of the data set (1,1,2,3,5). Left, center, and right are for $\bar r = (2, 1, .2)$.

Because the second derivative is not positive everywhere, the Cauchy function introduces the possibility of multiple solutions, but because of the good results we see in Figure 4, you might like to try it anyway. Perhaps the reason it seems to work so well is that it uses mostly residuals of ``average size,'' not the big ones or the small ones. This happens because $\Delta \bold m$ is made from $\bold F'$ and the components of $\bold W^2\bold r$ which are a function $r_i/(1+r_i^2/\bar r^2)$that is maximum for those residuals near $\bar r$.

Module irls supplies two useful weighting functions that can be interchanged as arguments to the reweighted scheme .  

#$=head1 NAME
#$
#$irls  - weighting functions
#$
#$=head1 SYNOPSIS
#$
#$C<ierr=l1(res,weight)>
#$
#$C<ierr=cauchy(res,weight)>
#$
#$=head1 PARAMETERS
#$
#$=over 4
#$
#$=item res  -  C<real(:)>  
#$
#$      Residual
#$
#$=item weight  -  C<real(:)>
#$
#$      Resulting weighting function
#$
#$=back
#$
#$=head1 DESCRIPTION
#$
#$Weighting function to apply to cgmethods (warning stability
#$is not guaranteed when using these).
#$
#$L1 -  weight=1/abs(1+ res/rbar)
#$CAUCHY- weight=1/sqrt(1+res/rbar**2)
#$
#$
#$=head1 SEE ALSO
#$
#$L<solver>
#$
#$=head1 LIBRARY
#$
#$B<geef90>
#$
#$=cut
module irls {
  use quantile_mod
contains
  integer function l1 (res, weight)  {
    real, dimension (:)  :: res, weight
    real                 :: rbar
    rbar = quantile( int( 0.5*size(res)), abs (res))            # median
    weight = 1. / sqrt( sqrt (1. + (res/rbar)**2));   l1 = 0
    }
  integer function cauchy (res, weight)  {
    real, dimension (:)  :: res, weight
    real                 :: rbar
    rbar = quantile( int( 0.5*size(res)), abs (res))             # median
    weight = 1. / sqrt (1. + (res/rbar)**2);          cauchy = 0
    }
}