Source Code for Module gmisclib.matrix_arrange

  1  """Take a covariance-like matrix, and re-order it to be close to a diagonal matrix. 
  2  We calculate a map that renames the indices (the same map for both left and right), 
  3  to bring large off-diagonal elements close to the diagonal. 
  4  """ 
  5   
  6  import math 
  7  import random 
  8  import numpy 
  9   
 10   
 11 -def _eval_map(a, m, wt): 
 12          """This is a measure of how far the matrix is from diagonal, 
 13          after application of the mapping. 
 14          Matrix "a" is a numeric python 2-d array. 
 15          Map is a numeric python 1-d integer array.""" 
 16   
 17          n = a.shape[0] 
 18          assert a.shape == (n, n) 
 19          # print "mapped a=", map_array(a, m) 
 20          # print "m=", m 
 21          return numpy.sum(map_array(a, m)*wt) 
 22   
 23   
 24 -def _offdiag(n): 
 25          index = numpy.arange(n) 
 26          offdiag = numpy.absolute(index[:,numpy.newaxis]-index[numpy.newaxis,:]) 
 27          assert offdiag.shape == (n, n) 
 28          assert offdiag[0, 0] == 0 
 29          assert offdiag[1, 0] == 1 
 30          # print "WT=", offdiag 
 31          return offdiag 
 32   
 33   
 34   
 35 -def _optimistic(n): 
 36          return numpy.ones((n, n), numpy.int) - numpy.identity(n) 
 37   
 38   
 39   
 40 -def diagonalize(a0): 
 41          EPS = 1e-10 
 42          a = numpy.array(a0, copy=True) 
 43          n = a.shape[0] 
 44          assert a.shape == (n, n) 
 45          wt = numpy.square(_offdiag(n)) 
 46          m = numpy.arange(n) 
 47          opval = _eval_map(a, m, _optimistic(n)) 
 48          # print "OPVAL=", opval 
 49          bad = _eval_map(a, m, wt) 
 50          bad0 = bad 
 51          f = 1.0 
 52          fr = math.pow(0.5, 1.0/n) 
 53          nobetter = 0 
 54          while 1: 
 55                  # i = random.choice(_worst(a, m, wt)) 
 56                  i = random.randrange(n) 
 57                  j = random.randrange(n) 
 58                  if j == i: 
 59                          continue 
 60                  # print "TESTING", i, j 
 61                  tmp = m[i] 
 62                  m[i] = m[j] 
 63                  m[j] = tmp 
 64                  bt = _eval_map(a, m, wt) 
 65                  if bt <= opval + bad*EPS: 
 66                          bad = bt 
 67                          break 
 68                  if bt <= bad * (1 + f/float(n)): 
 69                          # print "ACCEPT", bad, bt 
 70                          if abs(bt-bad) > 0.1*bad*f/float(n): 
 71                                  f *= fr 
 72                                  nobetter = 0 
 73                          bad = bt 
 74                          nobetter = 0 
 75                  else: 
 76                          # print "NG", bad, bt 
 77                          tmp = m[i] 
 78                          m[i] = m[j] 
 79                          m[j] = tmp 
 80                          nobetter += 1 
 81   
 82                  if nobetter >= n*n: 
 83                          break 
 84          if m[0] > m[-1]:        # Pick one of two equivalent orderings 
 85                  m = m[::-1] 
 86          return (m, bad0, bad) 
 87   
 88   
 89   
 90 -def map_array(a, m): 
 91          return numpy.take(numpy.take(a, m, axis=0), m, axis=1) 
 92   
 93   
 94 -def test(): 
 95          a = numpy.array([[1, 0, 3, 2], [0, 1, 0, 0], [3, 0, 1, 0], [2, 0, 0, 1]]) 
 96          m, b0, b1 = diagonalize(a) 
 97          assert b1<b0 
 98          print m 
 99          print map_array(a, m) 
100   
101   
102  if __name__ == '__main__': 
103          test() 
104