[SciPy-dev] Type handling of matrices

[Alexander Schmolck]

Nils Wagner <nwagner at mecha.uni-stuttgart.de> writes:

> BTW, how can I compute the gap function of A in scipy, where A is a dense
> normal matrix
> and gap(A) is defined as follows
>
> gap(A) = \min\limits_{i \ne j} | \lambda_i-\lambda_j | / 2
>
> This approach seems to be not very efficient
>
> def gap(A):
>  w = linalg.eigvals(A) # Compute the spectrum of A
>  min = 1.e10 # Initialize
>  for i in len(w):
>   for j in len(w):
>         if i <> j:
>               min1 = abs(w[i]-w[j])
>               if min1 < min:
>                   min=min1
>  return min/2

How about something along these lines (untested)?

def gap(A):
    w = linalg.eigvals(A)
    dt = abs(subtract.outer(w, w)).flat
    dt[::len(w)] = inf
    return minimum.reduce(dt) / 2.0
##    or altnernatively something like:
#    return minimum.reduce(nonzero( * 
#       not_equal.outer(arange(len(w)), arange(len(w))))).flat) / 2.0
    
'as