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, gesdd, 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
 144     # Create a copy of ``symmat`` here because ``syevr`` clobbers it.
 145     dummy = matrix(symmat, symmat.size)
 146     syevr(dummy, eigs)
 147     return list(eigs)
 148
 149
 150 def eigenvalues_re(anymat):
 151     """
 152     Return the real parts of the eigenvalues of the given square matrix.
 153
 154     Parameters
 155     ----------
 156
 157     anymat : matrix
 158         The square matrix whose eigenvalues you want.
 159
 160     Returns
 161     -------
 162
 163     list of float
 164        A list of the real parts (in no particular order) of the
 165        eigenvalues of ``anymat``.
 166
 167     Raises
 168     ------
 169
 170     TypeError
 171         If the input matrix is not square.
 172
 173     Examples
 174     --------
 175
 176     This is symmetric and has two real eigenvalues:
 177
 178         >>> A = matrix([[2,1],[1,2]], tc='d')
 179         >>> sorted(eigenvalues_re(A))
 180         [1.0, 3.0]
 181
 182     But this rotation matrix has eigenvalues `i` and `-i`, both of whose
 183     real parts are zero:
 184
 185         >>> A = matrix([[0,-1],[1,0]])
 186         >>> eigenvalues_re(A)
 187         [0.0, 0.0]
 188
 189     If the input matrix is not square, it doesn't have eigenvalues::
 190
 191         >>> A = matrix([[1,2],[3,4],[5,6]])
 192         >>> eigenvalues_re(A)
 193         Traceback (most recent call last):
 194         ...
 195         TypeError: input matrix must be square
 196
 197     """
 198     if not anymat.size[0] == anymat.size[1]:
 199         raise TypeError('input matrix must be square')
 200
 201     domain_dim = anymat.size[0]
 202     eigs = matrix(0, (domain_dim, 1), tc='z')
 203
 204     # Create a copy of ``anymat`` here for two reasons:
 205     #
 206     #   1. ``gees`` clobbers its input.
 207     #   2. We need to ensure that the type code of ``dummy`` is 'd' or 'z'.
 208     #
 209     dummy = matrix(anymat, anymat.size, tc='d')
 210
 211     gees(dummy, eigs)
 212     return [eig.real for eig in eigs]
 213
 214
 215 def identity(domain_dim, typecode='i'):
 216     """
 217     Create an identity matrix of the given dimensions.
 218
 219     Parameters
 220     ----------
 221
 222     domain_dim : int
 223         The dimension of the vector space on which the identity will act.
 224
 225     typecode : {'i', 'd', 'z'}, optional
 226         The type code for the returned matrix, defaults to 'i' for integers.
 227         Can also be 'd' for real double, or 'z' for complex double.
 228
 229     Returns
 230     -------
 231
 232     matrix
 233         A ``domain_dim``-by-``domain_dim`` dense integer identity matrix.
 234
 235     Raises
 236     ------
 237
 238     ValueError
 239         If you ask for the identity on zero or fewer dimensions.
 240
 241     Examples
 242     --------
 243
 244        >>> print(identity(3))
 245        [ 1  0  0]
 246        [ 0  1  0]
 247        [ 0  0  1]
 248        <BLANKLINE>
 249
 250     """
 251     if domain_dim <= 0:
 252         raise ValueError('domain dimension must be positive')
 253
 254     entries = [int(i == j)
 255                for i in range(domain_dim)
 256                for j in range(domain_dim)]
 257     return matrix(entries, (domain_dim, domain_dim), tc=typecode)
 258
 259
 260 def inner_product(vec1, vec2):
 261     """
 262     Compute the Euclidean inner product of two vectors.
 263
 264     Parameters
 265     ----------
 266
 267     vec1, vec2 : matrix
 268         The two vectors whose inner product you want.
 269
 270     Returns
 271     -------
 272
 273     float
 274         The inner product of ``vec1`` and ``vec2``.
 275
 276     Raises
 277     ------
 278
 279     TypeError
 280         If the lengths of ``vec1`` and ``vec2`` differ.
 281
 282     Examples
 283     --------
 284
 285         >>> x = [1,2,3]
 286         >>> y = [3,4,1]
 287         >>> inner_product(x,y)
 288         14
 289
 290         >>> x = matrix([1,1,1])
 291         >>> y = matrix([2,3,4], (1,3))
 292         >>> inner_product(x,y)
 293         9
 294
 295         >>> x = [1,2,3]
 296         >>> y = [1,1]
 297         >>> inner_product(x,y)
 298         Traceback (most recent call last):
 299         ...
 300         TypeError: the lengths of vec1 and vec2 must match
 301
 302     """
 303     if not len(vec1) == len(vec2):
 304         raise TypeError('the lengths of vec1 and vec2 must match')
 305
 306     return sum([x*y for (x, y) in zip(vec1, vec2)])
 307
 308
 309 def norm(matrix_or_vector):
 310     """
 311     Return the Frobenius norm of a matrix or vector.
 312
 313     When the input is a vector, its matrix-Frobenius norm is the same
 314     thing as its vector-Euclidean norm.
 315
 316     Parameters
 317     ----------
 318
 319     matrix_or_vector : matrix
 320         The matrix or vector whose norm you want.
 321
 322     Returns
 323     -------
 324
 325     float
 326         The norm of ``matrix_or_vector``.
 327
 328     Examples
 329     --------
 330
 331         >>> v = matrix([1,1])
 332         >>> print('{:.5f}'.format(norm(v)))
 333         1.41421
 334
 335         >>> A = matrix([1,1,1,1], (2,2))
 336         >>> norm(A)
 337         2.0
 338
 339     """
 340     return sqrt(inner_product(matrix_or_vector, matrix_or_vector))
 341
 342
 343 def vec(mat):
 344     """
 345     Create a long vector in column-major order from ``mat``.
 346
 347     Parameters
 348     ----------
 349
 350     mat : matrix
 351         Any sort of real matrix that you want written as a long vector.
 352
 353     Returns
 354     -------
 355
 356     matrix
 357         An ``len(mat)``-by-``1`` long column vector containign the
 358         entries of ``mat`` in column major order.
 359
 360     Examples
 361     --------
 362
 363         >>> A = matrix([[1,2],[3,4]])
 364         >>> print(A)
 365         [ 1  3]
 366         [ 2  4]
 367         <BLANKLINE>
 368
 369         >>> print(vec(A))
 370         [ 1]
 371         [ 2]
 372         [ 3]
 373         [ 4]
 374         <BLANKLINE>
 375
 376     Note that if ``mat`` is a vector, this function is a no-op:
 377
 378         >>> v = matrix([1,2,3,4], (4,1))
 379         >>> print(v)
 380         [ 1]
 381         [ 2]
 382         [ 3]
 383         [ 4]
 384         <BLANKLINE>
 385         >>> print(vec(v))
 386         [ 1]
 387         [ 2]
 388         [ 3]
 389         [ 4]
 390         <BLANKLINE>
 391
 392     """
 393     return matrix(mat, (len(mat), 1))
 394
 395
 396 def condition_number(mat):
 397     """
 398     Return the condition number of the given matrix.
 399
 400     The condition number of a matrix quantifies how hard it is to do
 401     numerical computation with that matrix. It is usually defined as
 402     the ratio of the norm of the matrix to the norm of its inverse, and
 403     therefore depends on the norm used. One way to compute the condition
 404     number with respect to the 2-norm is as the ratio of the matrix's
 405     largest and smallest singular values. Since we have easy access to
 406     those singular values, that is the algorithm we use.
 407
 408     The larger the condition number is, the worse the matrix is.
 409
 410     Parameters
 411     ----------
 412     mat : matrix
 413         The matrix whose condition number you want.
 414
 415     Returns
 416     -------
 417
 418     float
 419         The nonnegative condition number of ``mat``.
 420
 421     Examples
 422     --------
 423
 424     >>> condition_number(identity(1, typecode='d'))
 425     1.0
 426     >>> condition_number(identity(2, typecode='d'))
 427     1.0
 428     >>> condition_number(identity(3, typecode='d'))
 429     1.0
 430
 431     >>> A = matrix([[2,1],[1,2]], tc='d')
 432     >>> abs(condition_number(A) - 3.0) < options.ABS_TOL
 433     True
 434
 435     >>> A = matrix([[2,1j],[-1j,2]], tc='z')
 436     >>> abs(condition_number(A) - 3.0) < options.ABS_TOL
 437     True
 438
 439     """
 440     num_eigs = min(mat.size)
 441     eigs = matrix(0, (num_eigs, 1), tc='d')
 442     gesdd(mat, eigs)
 443
 444     if len(eigs) > 0:
 445         return eigs[0]/eigs[-1]
 446     else:
 447         return 0