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