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