[SciPy-user] Matrix sign function

[Pearu Peterson]

On Wed, 24 Sep 2003, Nils Wagner wrote:

> I have used linalg.signm with a defective matrix.  

The current algorithm in linalg.funm is not suitable for 
defective matrices. Btw, both Matlab and Octave funm functions
have the same problem. See for instance the corresponding thread in 
Octave mailing list:
  http://www.octave.org/octave-lists/archive/bug-octave.2001/msg00063.html

> The results do not fulfill the properties of a certain projector.  
> Please find enclosed my test program.  
>   Any comments or suggestions for a further improvement of signm ?

Higham in his Matlab Matrix Computation Toolbox uses similar algorithm to
linalg.funm (with certain steps changed for stability) so that I doubt 
that it would resolve the problem for defective matrices (correct me if I 
am wrong).

But using iterative method (from the references you gave), please test the 
following signm function with your matrices:

#---------------
feps = linalg.matfuncs.feps
eps = linalg.matfuncs.eps
_array_precision = linalg.matfuncs._array_precision
norm = linalg.norm
def signm3(a):
    """matrix sign"""
    def rounded_sign(x):
        rx = real(x)
        mx = amax(absolute(x)) # XXX: amax segfaults on complex input
        if rx.typecode()=='f':
            c =  1e3*feps*mx
        else:
            c =  1e3*eps*mx
        return sign( (absolute(rx) > c) * rx )
    result,err = funm(a, rounded_sign, disp=0)
    tol = {0:feps, 1:eps}[_array_precision[result.typecode()]]
    if err < 1000*tol:
        return result
    S0 = result + 0.5*identity(result.shape[0])/norm(result,1)
    for i in range(20):
        iS0 = linalg.inv(S0)
        alpha = beta = 0.5
        S0 = alpha*S0 + beta*iS0
        Pp=0.5*(dot(S0,S0)+S0)
        eps1 = norm(dot(Pp,Pp)-Pp)
        if eps1 < 1e3*tol:
            break
    return S0
#--------------

Pearu