src/dunshire/matrices.py

   1 """
   2 Utility functions for working with CVXOPT matrices (instances of the
   3 ``cvxopt.base.matrix`` class).
   4 """
   5
   6 from math import sqrt
   7 from cvxopt import matrix
   8 from cvxopt.lapack import syev
   9
  10 def append_col(left, right):
  11     """
  12     Append the matrix ``right`` to the right side of the matrix ``left``.
  13
  14     EXAMPLES:
  15
  16         >>> A = matrix([1,2,3,4], (2,2))
  17         >>> B = matrix([5,6,7,8,9,10], (2,3))
  18         >>> print(append_col(A,B))
  19         [  1   3   5   7   9]
  20         [  2   4   6   8  10]
  21         <BLANKLINE>
  22
  23     """
  24     return matrix([left.trans(), right.trans()]).trans()
  25
  26 def append_row(top, bottom):
  27     """
  28     Append the matrix ``bottom`` to the bottom of the matrix ``top``.
  29
  30     EXAMPLES:
  31
  32         >>> A = matrix([1,2,3,4], (2,2))
  33         >>> B = matrix([5,6,7,8,9,10], (3,2))
  34         >>> print(append_row(A,B))
  35         [  1   3]
  36         [  2   4]
  37         [  5   8]
  38         [  6   9]
  39         [  7  10]
  40         <BLANKLINE>
  41
  42     """
  43     return matrix([top, bottom])
  44
  45
  46 def eigenvalues(real_matrix):
  47     """
  48     Return the eigenvalues of the given ``real_matrix``.
  49
  50     EXAMPLES:
  51
  52         >>> A = matrix([[2,1],[1,2]], tc='d')
  53         >>> eigenvalues(A)
  54         [1.0, 3.0]
  55
  56     """
  57     domain_dim = real_matrix.size[0] # Assume ``real_matrix`` is square.
  58     eigs = matrix(0, (domain_dim, 1), tc='d')
  59     syev(real_matrix, eigs)
  60     return list(eigs)
  61
  62
  63 def identity(domain_dim):
  64     """
  65     Return a ``domain_dim``-by-``domain_dim`` dense integer identity
  66     matrix.
  67
  68     EXAMPLES:
  69
  70        >>> print(identity(3))
  71        [ 1  0  0]
  72        [ 0  1  0]
  73        [ 0  0  1]
  74        <BLANKLINE>
  75
  76     """
  77     if domain_dim <= 0:
  78         raise ValueError('domain dimension must be positive')
  79
  80     entries = [int(i == j)
  81                for i in range(domain_dim)
  82                for j in range(domain_dim)]
  83     return matrix(entries, (domain_dim, domain_dim))
  84
  85
  86 def inner_product(vec1, vec2):
  87     """
  88     Compute the (Euclidean) inner product of the two vectors ``vec1``
  89     and ``vec2``.
  90
  91     EXAMPLES:
  92
  93         >>> x = [1,2,3]
  94         >>> y = [3,4,1]
  95         >>> inner_product(x,y)
  96         14
  97
  98         >>> x = matrix([1,1,1])
  99         >>> y = matrix([2,3,4], (1,3))
 100         >>> inner_product(x,y)
 101         9
 102
 103         >>> x = [1,2,3]
 104         >>> y = [1,1]
 105         >>> inner_product(x,y)
 106         Traceback (most recent call last):
 107         ...
 108         TypeError: the lengths of vec1 and vec2 must match
 109
 110     """
 111     if not len(vec1) == len(vec2):
 112         raise TypeError('the lengths of vec1 and vec2 must match')
 113
 114     return sum([x*y for (x,y) in zip(vec1,vec2)])
 115
 116
 117 def norm(matrix_or_vector):
 118     """
 119     Return the Frobenius norm of ``matrix_or_vector``, which is the same
 120     thing as its Euclidean norm when it's a vector (when one of its
 121     dimensions is unity).
 122
 123     EXAMPLES:
 124
 125         >>> v = matrix([1,1])
 126         >>> print('{:.5f}'.format(norm(v)))
 127         1.41421
 128
 129         >>> A = matrix([1,1,1,1], (2,2))
 130         >>> norm(A)
 131         2.0
 132
 133     """
 134     return sqrt(inner_product(matrix_or_vector,matrix_or_vector))
 135
 136
 137 def vec(real_matrix):
 138     """
 139     Create a long vector in column-major order from ``real_matrix``.
 140
 141     EXAMPLES:
 142
 143         >>> A = matrix([[1,2],[3,4]])
 144         >>> print(A)
 145         [ 1  3]
 146         [ 2  4]
 147         <BLANKLINE>
 148
 149         >>> print(vec(A))
 150         [ 1]
 151         [ 2]
 152         [ 3]
 153         [ 4]
 154         <BLANKLINE>
 155
 156     Note that if ``real_matrix`` is a vector, this function is a no-op:
 157
 158         >>> v = matrix([1,2,3,4], (4,1))
 159         >>> print(v)
 160         [ 1]
 161         [ 2]
 162         [ 3]
 163         [ 4]
 164         <BLANKLINE>
 165         >>> print(vec(v))
 166         [ 1]
 167         [ 2]
 168         [ 3]
 169         [ 4]
 170         <BLANKLINE>
 171
 172     """
 173     return matrix(real_matrix, (len(real_matrix), 1))