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 gees, syevr
   9
  10 from . 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     On the surface, this appears redundant to the :func:`eigenvalues_re`
  99     function. However, if we know in advance that our input is
 100     symmetric, a better algorithm can be used.
 101
 102     Parameters
 103     ----------
 104
 105     symmat : matrix
 106         The real symmetric matrix whose eigenvalues you want.
 107
 108     Returns
 109     -------
 110
 111     list of float
 112        A list of the eigenvalues (in no particular order) of ``symmat``.
 113
 114     Raises
 115     ------
 116
 117     TypeError
 118         If the input matrix is not symmetric.
 119
 120     Examples
 121     --------
 122
 123         >>> A = matrix([[2,1],[1,2]], tc='d')
 124         >>> eigenvalues(A)
 125         [1.0, 3.0]
 126
 127     If the input matrix is not symmetric, it may not have real
 128     eigenvalues, and we don't know what to do::
 129
 130         >>> A = matrix([[1,2],[3,4]])
 131         >>> eigenvalues(A)
 132         Traceback (most recent call last):
 133         ...
 134         TypeError: input must be a symmetric real matrix
 135
 136     """
 137     if not norm(symmat.trans() - symmat) < options.ABS_TOL:
 138         # Ensure that ``symmat`` is symmetric (and thus square).
 139         raise TypeError('input must be a symmetric real matrix')
 140
 141     domain_dim = symmat.size[0]
 142     eigs = matrix(0, (domain_dim, 1), tc='d')
 143     syevr(symmat, eigs)
 144     return list(eigs)
 145
 146
 147 def eigenvalues_re(anymat):
 148     """
 149     Return the real parts of the eigenvalues of the given square matrix.
 150
 151     Parameters
 152     ----------
 153
 154     anymat : matrix
 155         The square matrix whose eigenvalues you want.
 156
 157     Returns
 158     -------
 159
 160     list of float
 161        A list of the real parts (in no particular order) of the
 162        eigenvalues of ``anymat``.
 163
 164     Raises
 165     ------
 166
 167     TypeError
 168         If the input matrix is not square.
 169
 170     Examples
 171     --------
 172
 173     This is symmetric and has two real eigenvalues:
 174
 175         >>> A = matrix([[2,1],[1,2]], tc='d')
 176         >>> sorted(eigenvalues_re(A))
 177         [1.0, 3.0]
 178
 179     But this rotation matrix has eigenvalues `i` and `-i`, both of whose
 180     real parts are zero:
 181
 182         >>> A = matrix([[0,-1],[1,0]])
 183         >>> eigenvalues_re(A)
 184         [0.0, 0.0]
 185
 186     If the input matrix is not square, it doesn't have eigenvalues::
 187
 188         >>> A = matrix([[1,2],[3,4],[5,6]])
 189         >>> eigenvalues_re(A)
 190         Traceback (most recent call last):
 191         ...
 192         TypeError: input matrix must be square
 193
 194     """
 195     if not anymat.size[0] == anymat.size[1]:
 196         raise TypeError('input matrix must be square')
 197
 198     domain_dim = anymat.size[0]
 199     eigs = matrix(0, (domain_dim, 1), tc='z')
 200
 201     # Create a copy of ``anymat`` here for two reasons:
 202     #
 203     #   1. ``gees`` clobbers its input.
 204     #   2. We need to ensure that the type code of ``dummy`` is 'd' or 'z'.
 205     #
 206     dummy = matrix(anymat, anymat.size, tc='d')
 207
 208     gees(dummy, eigs)
 209     return [eig.real for eig in eigs]
 210
 211
 212 def identity(domain_dim):
 213     """
 214     Create an identity matrix of the given dimensions.
 215
 216     Parameters
 217     ----------
 218
 219     domain_dim : int
 220         The dimension of the vector space on which the identity will act.
 221
 222     Returns
 223     -------
 224
 225     matrix
 226         A ``domain_dim``-by-``domain_dim`` dense integer identity matrix.
 227
 228     Raises
 229     ------
 230
 231     ValueError
 232         If you ask for the identity on zero or fewer dimensions.
 233
 234     Examples
 235     --------
 236
 237        >>> print(identity(3))
 238        [ 1  0  0]
 239        [ 0  1  0]
 240        [ 0  0  1]
 241        <BLANKLINE>
 242
 243     """
 244     if domain_dim <= 0:
 245         raise ValueError('domain dimension must be positive')
 246
 247     entries = [int(i == j)
 248                for i in range(domain_dim)
 249                for j in range(domain_dim)]
 250     return matrix(entries, (domain_dim, domain_dim))
 251
 252
 253 def inner_product(vec1, vec2):
 254     """
 255     Compute the Euclidean inner product of two vectors.
 256
 257     Parameters
 258     ----------
 259
 260     vec1, vec2 : matrix
 261         The two vectors whose inner product you want.
 262
 263     Returns
 264     -------
 265
 266     float
 267         The inner product of ``vec1`` and ``vec2``.
 268
 269     Raises
 270     ------
 271
 272     TypeError
 273         If the lengths of ``vec1`` and ``vec2`` differ.
 274
 275     Examples
 276     --------
 277
 278         >>> x = [1,2,3]
 279         >>> y = [3,4,1]
 280         >>> inner_product(x,y)
 281         14
 282
 283         >>> x = matrix([1,1,1])
 284         >>> y = matrix([2,3,4], (1,3))
 285         >>> inner_product(x,y)
 286         9
 287
 288         >>> x = [1,2,3]
 289         >>> y = [1,1]
 290         >>> inner_product(x,y)
 291         Traceback (most recent call last):
 292         ...
 293         TypeError: the lengths of vec1 and vec2 must match
 294
 295     """
 296     if not len(vec1) == len(vec2):
 297         raise TypeError('the lengths of vec1 and vec2 must match')
 298
 299     return sum([x*y for (x, y) in zip(vec1, vec2)])
 300
 301
 302 def norm(matrix_or_vector):
 303     """
 304     Return the Frobenius norm of a matrix or vector.
 305
 306     When the input is a vector, its matrix-Frobenius norm is the same
 307     thing as its vector-Euclidean norm.
 308
 309     Parameters
 310     ----------
 311
 312     matrix_or_vector : matrix
 313         The matrix or vector whose norm you want.
 314
 315     Returns
 316     -------
 317
 318     float
 319         The norm of ``matrix_or_vector``.
 320
 321     Examples
 322     --------
 323
 324         >>> v = matrix([1,1])
 325         >>> print('{:.5f}'.format(norm(v)))
 326         1.41421
 327
 328         >>> A = matrix([1,1,1,1], (2,2))
 329         >>> norm(A)
 330         2.0
 331
 332     """
 333     return sqrt(inner_product(matrix_or_vector, matrix_or_vector))
 334
 335
 336 def vec(mat):
 337     """
 338     Create a long vector in column-major order from ``mat``.
 339
 340     Parameters
 341     ----------
 342
 343     mat : matrix
 344         Any sort of real matrix that you want written as a long vector.
 345
 346     Returns
 347     -------
 348
 349     matrix
 350         An ``len(mat)``-by-``1`` long column vector containign the
 351         entries of ``mat`` in column major order.
 352
 353     Examples
 354     --------
 355
 356         >>> A = matrix([[1,2],[3,4]])
 357         >>> print(A)
 358         [ 1  3]
 359         [ 2  4]
 360         <BLANKLINE>
 361
 362         >>> print(vec(A))
 363         [ 1]
 364         [ 2]
 365         [ 3]
 366         [ 4]
 367         <BLANKLINE>
 368
 369     Note that if ``mat`` is a vector, this function is a no-op:
 370
 371         >>> v = matrix([1,2,3,4], (4,1))
 372         >>> print(v)
 373         [ 1]
 374         [ 2]
 375         [ 3]
 376         [ 4]
 377         <BLANKLINE>
 378         >>> print(vec(v))
 379         [ 1]
 380         [ 2]
 381         [ 3]
 382         [ 4]
 383         <BLANKLINE>
 384
 385     """
 386     return matrix(mat, (len(mat), 1))