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 norm(matrix_or_vector):
  87     """
  88     Return the Frobenius norm of ``matrix_or_vector``, which is the same
  89     thing as its Euclidean norm when it's a vector (when one of its
  90     dimensions is unity).
  91
  92     EXAMPLES:
  93
  94         >>> v = matrix([1,1])
  95         >>> print('{:.5f}'.format(norm(v)))
  96         1.41421
  97
  98         >>> A = matrix([1,1,1,1], (2,2))
  99         >>> norm(A)
 100         2.0
 101
 102     """
 103     return sqrt(sum([x**2 for x in matrix_or_vector]))
 104
 105
 106 def vec(real_matrix):
 107     """
 108     Create a long vector in column-major order from ``real_matrix``.
 109
 110     EXAMPLES:
 111
 112         >>> A = matrix([[1,2],[3,4]])
 113         >>> print(A)
 114         [ 1  3]
 115         [ 2  4]
 116         <BLANKLINE>
 117
 118         >>> print(vec(A))
 119         [ 1]
 120         [ 2]
 121         [ 3]
 122         [ 4]
 123         <BLANKLINE>
 124
 125     Note that if ``real_matrix`` is a vector, this function is a no-op:
 126
 127         >>> v = matrix([1,2,3,4], (4,1))
 128         >>> print(v)
 129         [ 1]
 130         [ 2]
 131         [ 3]
 132         [ 4]
 133         <BLANKLINE>
 134         >>> print(vec(v))
 135         [ 1]
 136         [ 2]
 137         [ 3]
 138         [ 4]
 139         <BLANKLINE>
 140
 141     """
 142     return matrix(real_matrix, (len(real_matrix), 1))