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     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, typecode='i'):
 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     typecode : {'i', 'd', 'z'}, optional
 223         The type code for the returned matrix, defaults to 'i' for integers.
 224         Can also be 'd' for real double, or 'z' for complex double.
 225
 226     Returns
 227     -------
 228
 229     matrix
 230         A ``domain_dim``-by-``domain_dim`` dense integer identity matrix.
 231
 232     Raises
 233     ------
 234
 235     ValueError
 236         If you ask for the identity on zero or fewer dimensions.
 237
 238     Examples
 239     --------
 240
 241        >>> print(identity(3))
 242        [ 1  0  0]
 243        [ 0  1  0]
 244        [ 0  0  1]
 245        <BLANKLINE>
 246
 247     """
 248     if domain_dim <= 0:
 249         raise ValueError('domain dimension must be positive')
 250
 251     entries = [int(i == j)
 252                for i in range(domain_dim)
 253                for j in range(domain_dim)]
 254     return matrix(entries, (domain_dim, domain_dim), tc=typecode)
 255
 256
 257 def inner_product(vec1, vec2):
 258     """
 259     Compute the Euclidean inner product of two vectors.
 260
 261     Parameters
 262     ----------
 263
 264     vec1, vec2 : matrix
 265         The two vectors whose inner product you want.
 266
 267     Returns
 268     -------
 269
 270     float
 271         The inner product of ``vec1`` and ``vec2``.
 272
 273     Raises
 274     ------
 275
 276     TypeError
 277         If the lengths of ``vec1`` and ``vec2`` differ.
 278
 279     Examples
 280     --------
 281
 282         >>> x = [1,2,3]
 283         >>> y = [3,4,1]
 284         >>> inner_product(x,y)
 285         14
 286
 287         >>> x = matrix([1,1,1])
 288         >>> y = matrix([2,3,4], (1,3))
 289         >>> inner_product(x,y)
 290         9
 291
 292         >>> x = [1,2,3]
 293         >>> y = [1,1]
 294         >>> inner_product(x,y)
 295         Traceback (most recent call last):
 296         ...
 297         TypeError: the lengths of vec1 and vec2 must match
 298
 299     """
 300     if not len(vec1) == len(vec2):
 301         raise TypeError('the lengths of vec1 and vec2 must match')
 302
 303     return sum([x*y for (x, y) in zip(vec1, vec2)])
 304
 305
 306 def norm(matrix_or_vector):
 307     """
 308     Return the Frobenius norm of a matrix or vector.
 309
 310     When the input is a vector, its matrix-Frobenius norm is the same
 311     thing as its vector-Euclidean norm.
 312
 313     Parameters
 314     ----------
 315
 316     matrix_or_vector : matrix
 317         The matrix or vector whose norm you want.
 318
 319     Returns
 320     -------
 321
 322     float
 323         The norm of ``matrix_or_vector``.
 324
 325     Examples
 326     --------
 327
 328         >>> v = matrix([1,1])
 329         >>> print('{:.5f}'.format(norm(v)))
 330         1.41421
 331
 332         >>> A = matrix([1,1,1,1], (2,2))
 333         >>> norm(A)
 334         2.0
 335
 336     """
 337     return sqrt(inner_product(matrix_or_vector, matrix_or_vector))
 338
 339
 340 def vec(mat):
 341     """
 342     Create a long vector in column-major order from ``mat``.
 343
 344     Parameters
 345     ----------
 346
 347     mat : matrix
 348         Any sort of real matrix that you want written as a long vector.
 349
 350     Returns
 351     -------
 352
 353     matrix
 354         An ``len(mat)``-by-``1`` long column vector containign the
 355         entries of ``mat`` in column major order.
 356
 357     Examples
 358     --------
 359
 360         >>> A = matrix([[1,2],[3,4]])
 361         >>> print(A)
 362         [ 1  3]
 363         [ 2  4]
 364         <BLANKLINE>
 365
 366         >>> print(vec(A))
 367         [ 1]
 368         [ 2]
 369         [ 3]
 370         [ 4]
 371         <BLANKLINE>
 372
 373     Note that if ``mat`` is a vector, this function is a no-op:
 374
 375         >>> v = matrix([1,2,3,4], (4,1))
 376         >>> print(v)
 377         [ 1]
 378         [ 2]
 379         [ 3]
 380         [ 4]
 381         <BLANKLINE>
 382         >>> print(vec(v))
 383         [ 1]
 384         [ 2]
 385         [ 3]
 386         [ 4]
 387         <BLANKLINE>
 388
 389     """
 390     return matrix(mat, (len(mat), 1))
 391
 392
 393 def condition_number(mat):
 394     """
 395     Return the condition number of the given matrix.
 396
 397     The condition number of a matrix quantifies how hard it is to do
 398     numerical computation with that matrix. It is usually defined as
 399     the ratio of the norm of the matrix to the norm of its inverse, and
 400     therefore depends on the norm used. One way to compute the condition
 401     number with respect to the 2-norm is as the ratio of the matrix's
 402     largest and smallest singular values. Since we have easy access to
 403     those singular values, that is the algorithm we use.
 404
 405     The larger the condition number is, the worse the matrix is.
 406
 407     Parameters
 408     ----------
 409     mat : matrix
 410         The matrix whose condition number you want.
 411
 412     Returns
 413     -------
 414
 415     float
 416         The nonnegative condition number of ``mat``.
 417
 418     Examples
 419     --------
 420
 421     >>> condition_number(identity(1, typecode='d'))
 422     1.0
 423     >>> condition_number(identity(2, typecode='d'))
 424     1.0
 425     >>> condition_number(identity(3, typecode='d'))
 426     1.0
 427
 428     >>> A = matrix([[2,1],[1,2]], tc='d')
 429     >>> abs(condition_number(A) - 3.0) < options.ABS_TOL
 430     True
 431
 432     >>> A = matrix([[2,1j],[-1j,2]], tc='z')
 433     >>> abs(condition_number(A) - 3.0) < options.ABS_TOL
 434     True
 435
 436     """
 437     num_eigs = min(mat.size)
 438     eigs = matrix(0, (num_eigs, 1), tc='d')
 439     gesdd(mat, eigs)
 440
 441     if len(eigs) > 0:
 442         return eigs[0]/eigs[-1]
 443     else:
 444         return 0