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