src/dunshire/matrices.py

   1 """
   2 Utility functions for working with CVXOPT matrices (instances of the
   3 class:`cvxopt.base.matrix` class).
   4 """
   5
   6 from math import sqrt
   7 from cvxopt import matrix
   8 from cvxopt.lapack import syev
   9
  10 import options
  11
  12 def append_col(left, right):
  13     """
  14     Append the matrix ``right`` to the right side of the matrix ``left``.
  15
  16     Parameters
  17     ----------
  18     left, right : matrix
  19         The two matrices to append to one another.
  20
  21     Examples
  22     --------
  23
  24         >>> A = matrix([1,2,3,4], (2,2))
  25         >>> B = matrix([5,6,7,8,9,10], (2,3))
  26         >>> print(A)
  27         [ 1  3]
  28         [ 2  4]
  29         <BLANKLINE>
  30         >>> print(B)
  31         [  5   7   9]
  32         [  6   8  10]
  33         <BLANKLINE>
  34         >>> print(append_col(A,B))
  35         [  1   3   5   7   9]
  36         [  2   4   6   8  10]
  37         <BLANKLINE>
  38
  39     """
  40     return matrix([left.trans(), right.trans()]).trans()
  41
  42 def append_row(top, bottom):
  43     """
  44     Append the matrix ``bottom`` to the bottom of the matrix ``top``.
  45
  46     Parameters
  47     ----------
  48     top, bottom : matrix
  49         The two matrices to append to one another.
  50
  51     Examples
  52     --------
  53
  54         >>> A = matrix([1,2,3,4], (2,2))
  55         >>> B = matrix([5,6,7,8,9,10], (3,2))
  56         >>> print(A)
  57         [ 1  3]
  58         [ 2  4]
  59         <BLANKLINE>
  60         >>> print(B)
  61         [  5   8]
  62         [  6   9]
  63         [  7  10]
  64         <BLANKLINE>
  65         >>> print(append_row(A,B))
  66         [  1   3]
  67         [  2   4]
  68         [  5   8]
  69         [  6   9]
  70         [  7  10]
  71         <BLANKLINE>
  72
  73     """
  74     return matrix([top, bottom])
  75
  76
  77 def eigenvalues(symmat):
  78     """
  79     Return the eigenvalues of the given symmetric real matrix.
  80
  81     Parameters
  82     ----------
  83
  84     symmat : matrix
  85         The real symmetric matrix whose eigenvalues you want.
  86
  87
  88     Raises
  89     ------
  90     TypeError
  91         If the input matrix is not symmetric.
  92
  93     Examples
  94     --------
  95
  96         >>> A = matrix([[2,1],[1,2]], tc='d')
  97         >>> eigenvalues(A)
  98         [1.0, 3.0]
  99
 100     If the input matrix is not symmetric, it may not have real
 101     eigenvalues, and we don't know what to do::
 102
 103         >>> A = matrix([[1,2],[3,4]])
 104         >>> eigenvalues(A)
 105         Traceback (most recent call last):
 106         ...
 107         TypeError: input must be a symmetric real matrix
 108
 109     """
 110     if not norm(symmat.trans() - symmat) < options.ABS_TOL:
 111         # Ensure that ``symmat`` is symmetric (and thus square).
 112         raise TypeError('input must be a symmetric real matrix')
 113
 114     domain_dim = symmat.size[0]
 115     eigs = matrix(0, (domain_dim, 1), tc='d')
 116     syev(symmat, eigs)
 117     return list(eigs)
 118
 119
 120 def identity(domain_dim):
 121     """
 122     Create an identity matrix of the given dimensions.
 123
 124     Parameters
 125     ----------
 126
 127     domain_dim : int
 128         The dimension of the vector space on which the identity will act.
 129
 130     Returns
 131     -------
 132     matrix
 133         A ``domain_dim``-by-``domain_dim`` dense integer identity matrix.
 134
 135     Examples
 136     --------
 137
 138        >>> print(identity(3))
 139        [ 1  0  0]
 140        [ 0  1  0]
 141        [ 0  0  1]
 142        <BLANKLINE>
 143
 144     """
 145     if domain_dim <= 0:
 146         raise ValueError('domain dimension must be positive')
 147
 148     entries = [int(i == j)
 149                for i in range(domain_dim)
 150                for j in range(domain_dim)]
 151     return matrix(entries, (domain_dim, domain_dim))
 152
 153
 154 def inner_product(vec1, vec2):
 155     """
 156     Compute the Euclidean inner product of two vectors.
 157
 158     Parameters
 159     ----------
 160
 161     vec1, vec2 : matrix
 162         The two vectors whose inner product you want.
 163
 164     Returns
 165     -------
 166     float
 167         The inner product of ``vec1`` and ``vec2``.
 168
 169     Examples
 170     --------
 171
 172         >>> x = [1,2,3]
 173         >>> y = [3,4,1]
 174         >>> inner_product(x,y)
 175         14
 176
 177         >>> x = matrix([1,1,1])
 178         >>> y = matrix([2,3,4], (1,3))
 179         >>> inner_product(x,y)
 180         9
 181
 182         >>> x = [1,2,3]
 183         >>> y = [1,1]
 184         >>> inner_product(x,y)
 185         Traceback (most recent call last):
 186         ...
 187         TypeError: the lengths of vec1 and vec2 must match
 188
 189     """
 190     if not len(vec1) == len(vec2):
 191         raise TypeError('the lengths of vec1 and vec2 must match')
 192
 193     return sum([x*y for (x, y) in zip(vec1, vec2)])
 194
 195
 196 def norm(matrix_or_vector):
 197     """
 198     Return the Frobenius norm of a matrix or vector.
 199
 200     When the input is a vector, its matrix-Frobenius norm is the same
 201     thing as its vector-Euclidean norm.
 202
 203     Parameters
 204     ----------
 205
 206     matrix_or_vector : matrix
 207         The matrix or vector whose norm you want.
 208
 209     Returns
 210     -------
 211
 212     float
 213         The norm of ``matrix_or_vector``.
 214
 215     Examples
 216     --------
 217
 218         >>> v = matrix([1,1])
 219         >>> print('{:.5f}'.format(norm(v)))
 220         1.41421
 221
 222         >>> A = matrix([1,1,1,1], (2,2))
 223         >>> norm(A)
 224         2.0
 225
 226     """
 227     return sqrt(inner_product(matrix_or_vector, matrix_or_vector))
 228
 229
 230 def vec(mat):
 231     """
 232     Create a long vector in column-major order from ``mat``.
 233
 234     Parameters
 235     ----------
 236
 237     mat : matrix
 238         Any sort of real matrix that you want written as a long vector.
 239
 240     Returns
 241     -------
 242
 243     matrix
 244         An ``len(mat)``-by-``1`` long column vector containign the
 245         entries of ``mat`` in column major order.
 246
 247     Examples
 248     --------
 249
 250         >>> A = matrix([[1,2],[3,4]])
 251         >>> print(A)
 252         [ 1  3]
 253         [ 2  4]
 254         <BLANKLINE>
 255
 256         >>> print(vec(A))
 257         [ 1]
 258         [ 2]
 259         [ 3]
 260         [ 4]
 261         <BLANKLINE>
 262
 263     Note that if ``mat`` is a vector, this function is a no-op:
 264
 265         >>> v = matrix([1,2,3,4], (4,1))
 266         >>> print(v)
 267         [ 1]
 268         [ 2]
 269         [ 3]
 270         [ 4]
 271         <BLANKLINE>
 272         >>> print(vec(v))
 273         [ 1]
 274         [ 2]
 275         [ 3]
 276         [ 4]
 277         <BLANKLINE>
 278
 279     """
 280     return matrix(mat, (len(mat), 1))