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 : matrix
  22         The "original" matrix, the one that will wind up on the left.
  23
  24     right : matrix
  25         The matrix to be appended on the right of ``left``.
  26
  27     Returns
  28     -------
  29
  30     matrix
  31         A new matrix consisting of ``right`` appended to the right
  32         of ``left``.
  33
  34     Examples
  35     --------
  36
  37         >>> A = matrix([1,2,3,4], (2,2))
  38         >>> B = matrix([5,6,7,8,9,10], (2,3))
  39         >>> print(A)
  40         [ 1  3]
  41         [ 2  4]
  42         <BLANKLINE>
  43         >>> print(B)
  44         [  5   7   9]
  45         [  6   8  10]
  46         <BLANKLINE>
  47         >>> print(append_col(A,B))
  48         [  1   3   5   7   9]
  49         [  2   4   6   8  10]
  50         <BLANKLINE>
  51
  52     """
  53     return matrix([left.trans(), right.trans()]).trans()
  54
  55
  56 def append_row(top, bottom):
  57     """
  58     Append two matrices top-to-bottom.
  59
  60     Parameters
  61     ----------
  62
  63     top : matrix
  64         The "original" matrix, the one that will wind up on top.
  65
  66     bottom : matrix
  67         The matrix to be appended below ``top``.
  68
  69     Returns
  70     -------
  71
  72     matrix
  73         A new matrix consisting of ``bottom`` appended below ``top``.
  74
  75     Examples
  76     --------
  77
  78         >>> A = matrix([1,2,3,4], (2,2))
  79         >>> B = matrix([5,6,7,8,9,10], (3,2))
  80         >>> print(A)
  81         [ 1  3]
  82         [ 2  4]
  83         <BLANKLINE>
  84         >>> print(B)
  85         [  5   8]
  86         [  6   9]
  87         [  7  10]
  88         <BLANKLINE>
  89         >>> print(append_row(A,B))
  90         [  1   3]
  91         [  2   4]
  92         [  5   8]
  93         [  6   9]
  94         [  7  10]
  95         <BLANKLINE>
  96
  97     """
  98     return matrix([top, bottom])
  99
 100
 101 def eigenvalues(symmat):
 102     """
 103     Return the eigenvalues of the given symmetric real matrix.
 104
 105     On the surface, this appears redundant to the :func:`eigenvalues_re`
 106     function. However, if we know in advance that our input is
 107     symmetric, a better algorithm can be used.
 108
 109     Parameters
 110     ----------
 111
 112     symmat : matrix
 113         The real symmetric matrix whose eigenvalues you want.
 114
 115     Returns
 116     -------
 117
 118     list of float
 119        A list of the eigenvalues (in no particular order) of ``symmat``.
 120
 121     Raises
 122     ------
 123
 124     TypeError
 125         If the input matrix is not symmetric.
 126
 127     Examples
 128     --------
 129
 130         >>> A = matrix([[2,1],[1,2]], tc='d')
 131         >>> eigenvalues(A)
 132         [1.0, 3.0]
 133
 134     If the input matrix is not symmetric, it may not have real
 135     eigenvalues, and we don't know what to do::
 136
 137         >>> A = matrix([[1,2],[3,4]])
 138         >>> eigenvalues(A)
 139         Traceback (most recent call last):
 140         ...
 141         TypeError: input must be a symmetric real matrix
 142
 143     """
 144     if not norm(symmat.trans() - symmat) < options.ABS_TOL:
 145         # Ensure that ``symmat`` is symmetric (and thus square).
 146         raise TypeError('input must be a symmetric real matrix')
 147
 148     domain_dim = symmat.size[0]
 149     eigs = matrix(0, (domain_dim, 1), tc='d')
 150
 151     # Create a copy of ``symmat`` here because ``syevr`` clobbers it.
 152     dummy = copy(symmat)
 153     syevr(dummy, eigs)
 154     return list(eigs)
 155
 156
 157 def eigenvalues_re(anymat):
 158     """
 159     Return the real parts of the eigenvalues of the given square matrix.
 160
 161     Parameters
 162     ----------
 163
 164     anymat : matrix
 165         The square matrix whose eigenvalues you want.
 166
 167     Returns
 168     -------
 169
 170     list of float
 171        A list of the real parts (in no particular order) of the
 172        eigenvalues of ``anymat``.
 173
 174     Raises
 175     ------
 176
 177     TypeError
 178         If the input matrix is not square.
 179
 180     Examples
 181     --------
 182
 183     This is symmetric and has two real eigenvalues:
 184
 185         >>> A = matrix([[2,1],[1,2]], tc='d')
 186         >>> sorted(eigenvalues_re(A))
 187         [1.0, 3.0]
 188
 189     But this rotation matrix has eigenvalues `i` and `-i`, both of whose
 190     real parts are zero:
 191
 192         >>> A = matrix([[0,-1],[1,0]])
 193         >>> eigenvalues_re(A)
 194         [0.0, 0.0]
 195
 196     If the input matrix is not square, it doesn't have eigenvalues::
 197
 198         >>> A = matrix([[1,2],[3,4],[5,6]])
 199         >>> eigenvalues_re(A)
 200         Traceback (most recent call last):
 201         ...
 202         TypeError: input matrix must be square
 203
 204     """
 205     if not anymat.size[0] == anymat.size[1]:
 206         raise TypeError('input matrix must be square')
 207
 208     domain_dim = anymat.size[0]
 209     eigs = matrix(0, (domain_dim, 1), tc='z')
 210
 211     # Create a copy of ``anymat`` here for two reasons:
 212     #
 213     #   1. ``gees`` clobbers its input.
 214     #   2. We need to ensure that the type code of ``dummy`` is 'd' or 'z'.
 215     #
 216     dummy = matrix(anymat, anymat.size, tc='d')
 217
 218     gees(dummy, eigs)
 219     return [eig.real for eig in eigs]
 220
 221
 222 def identity(domain_dim, typecode='i'):
 223     """
 224     Create an identity matrix of the given dimensions.
 225
 226     Parameters
 227     ----------
 228
 229     domain_dim : int
 230         The dimension of the vector space on which the identity will act.
 231
 232     typecode : {'i', 'd', 'z'}, optional
 233         The type code for the returned matrix, defaults to 'i' for integers.
 234         Can also be 'd' for real double, or 'z' for complex double.
 235
 236     Returns
 237     -------
 238
 239     matrix
 240         A ``domain_dim``-by-``domain_dim`` dense integer identity matrix.
 241
 242     Raises
 243     ------
 244
 245     ValueError
 246         If you ask for the identity on zero or fewer dimensions.
 247
 248     Examples
 249     --------
 250
 251        >>> print(identity(3))
 252        [ 1  0  0]
 253        [ 0  1  0]
 254        [ 0  0  1]
 255        <BLANKLINE>
 256
 257     """
 258     if domain_dim <= 0:
 259         raise ValueError('domain dimension must be positive')
 260
 261     entries = [int(i == j)
 262                for i in range(domain_dim)
 263                for j in range(domain_dim)]
 264     return matrix(entries, (domain_dim, domain_dim), tc=typecode)
 265
 266
 267 def inner_product(vec1, vec2):
 268     """
 269     Compute the Euclidean inner product of two vectors.
 270
 271     Parameters
 272     ----------
 273
 274     vec1 : matrix
 275         The first vector, whose inner product with ``vec2`` you want.
 276
 277     vec2 : matrix
 278         The second vector, whose inner product with ``vec1`` you want.
 279
 280     Returns
 281     -------
 282
 283     float
 284         The inner product of ``vec1`` and ``vec2``.
 285
 286     Raises
 287     ------
 288
 289     TypeError
 290         If the lengths of ``vec1`` and ``vec2`` differ.
 291
 292     Examples
 293     --------
 294
 295         >>> x = [1,2,3]
 296         >>> y = [3,4,1]
 297         >>> inner_product(x,y)
 298         14
 299
 300         >>> x = matrix([1,1,1])
 301         >>> y = matrix([2,3,4], (1,3))
 302         >>> inner_product(x,y)
 303         9
 304
 305         >>> x = [1,2,3]
 306         >>> y = [1,1]
 307         >>> inner_product(x,y)
 308         Traceback (most recent call last):
 309         ...
 310         TypeError: the lengths of vec1 and vec2 must match
 311
 312     """
 313     if not len(vec1) == len(vec2):
 314         raise TypeError('the lengths of vec1 and vec2 must match')
 315
 316     return sum([x*y for (x, y) in zip(vec1, vec2)])
 317
 318
 319 def norm(matrix_or_vector):
 320     """
 321     Return the Frobenius norm of a matrix or vector.
 322
 323     When the input is a vector, its matrix-Frobenius norm is the same
 324     thing as its vector-Euclidean norm.
 325
 326     Parameters
 327     ----------
 328
 329     matrix_or_vector : matrix
 330         The matrix or vector whose norm you want.
 331
 332     Returns
 333     -------
 334
 335     float
 336         The norm of ``matrix_or_vector``.
 337
 338     Examples
 339     --------
 340
 341         >>> v = matrix([1,1])
 342         >>> norm(v)
 343         1.414...
 344
 345         >>> A = matrix([1,1,1,1], (2,2))
 346         >>> norm(A)
 347         2.0...
 348
 349     """
 350     return sqrt(inner_product(matrix_or_vector, matrix_or_vector))
 351
 352
 353 def specnorm(mat):
 354     """
 355     Return the spectral norm of a matrix.
 356
 357     The spectral norm of a matrix is its largest singular value, and it
 358     corresponds to the operator norm induced by the vector Euclidean norm.
 359
 360     Parameters
 361     ----------
 362
 363     mat : matrix
 364         The matrix whose spectral norm you want.
 365
 366     Examples:
 367
 368         >>> specnorm(identity(3))
 369         1.0
 370
 371         >>> specnorm(5*identity(4))
 372         5.0
 373
 374     """
 375     num_eigs = min(mat.size)
 376     eigs = matrix(0, (num_eigs, 1), tc='d')
 377     typecode = 'd'
 378     if any([isinstance(entry, complex) for entry in mat]):
 379         typecode = 'z'
 380     dummy = matrix(mat, mat.size, tc=typecode)
 381     gesdd(dummy, eigs)
 382
 383     if len(eigs) > 0:
 384         return eigs[0]
 385     else:
 386         return 0
 387
 388
 389 def vec(mat):
 390     """
 391     Create a long vector in column-major order from ``mat``.
 392
 393     Parameters
 394     ----------
 395
 396     mat : matrix
 397         Any sort of real matrix that you want written as a long vector.
 398
 399     Returns
 400     -------
 401
 402     matrix
 403         An ``len(mat)``-by-``1`` long column vector containign the
 404         entries of ``mat`` in column major order.
 405
 406     Examples
 407     --------
 408
 409         >>> A = matrix([[1,2],[3,4]])
 410         >>> print(A)
 411         [ 1  3]
 412         [ 2  4]
 413         <BLANKLINE>
 414
 415         >>> print(vec(A))
 416         [ 1]
 417         [ 2]
 418         [ 3]
 419         [ 4]
 420         <BLANKLINE>
 421
 422     Note that if ``mat`` is a vector, this function is a no-op:
 423
 424         >>> v = matrix([1,2,3,4], (4,1))
 425         >>> print(v)
 426         [ 1]
 427         [ 2]
 428         [ 3]
 429         [ 4]
 430         <BLANKLINE>
 431         >>> print(vec(v))
 432         [ 1]
 433         [ 2]
 434         [ 3]
 435         [ 4]
 436         <BLANKLINE>
 437
 438     """
 439     return matrix(mat, (len(mat), 1))
 440
 441
 442 def condition_number(mat):
 443     """
 444     Return the condition number of the given matrix.
 445
 446     The condition number of a matrix quantifies how hard it is to do
 447     numerical computation with that matrix. It is usually defined as
 448     the ratio of the norm of the matrix to the norm of its inverse, and
 449     therefore depends on the norm used. One way to compute the condition
 450     number with respect to the 2-norm is as the ratio of the matrix's
 451     largest and smallest singular values. Since we have easy access to
 452     those singular values, that is the algorithm we use.
 453
 454     The larger the condition number is, the worse the matrix is.
 455
 456     Parameters
 457     ----------
 458     mat : matrix
 459         The matrix whose condition number you want.
 460
 461     Returns
 462     -------
 463
 464     float
 465         The nonnegative condition number of ``mat``.
 466
 467     Examples
 468     --------
 469
 470     >>> condition_number(identity(3))
 471     1.0
 472
 473     >>> A = matrix([[2,1],[1,2]])
 474     >>> abs(condition_number(A) - 3.0) < options.ABS_TOL
 475     True
 476
 477     >>> A = matrix([[2,1j],[-1j,2]])
 478     >>> abs(condition_number(A) - 3.0) < options.ABS_TOL
 479     True
 480
 481     """
 482     num_eigs = min(mat.size)
 483     eigs = matrix(0, (num_eigs, 1), tc='d')
 484     typecode = 'd'
 485     if any([isinstance(entry, complex) for entry in mat]):
 486         typecode = 'z'
 487     dummy = matrix(mat, mat.size, tc=typecode)
 488     gesdd(dummy, eigs)
 489
 490     if len(eigs) > 0:
 491         return eigs[0]/eigs[-1]
 492     else:
 493         return 0