[SciPy-user] Matrix sign function

[Nils Wagner]

------------------- 
 
 
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 
 
So, Scipy has the chance to be better than other tools. 
 
> 
> 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: 
 
Thank you very much for your modified signm function. 
I have used your signm function with the defective matrix. It works 
fine ! However, this a not a proof in general. 
 
Nils 
 
#--------------- 
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 
 
_______________________________________________ 
SciPy-user mailing list 
SciPy-user at scipy.net 
http://www.scipy.net/mailman/listinfo/scipy-user