]>
gitweb.michael.orlitzky.com - dunshire.git/blob - dunshire/matrices.py
2 Utility functions for working with CVXOPT matrices (instances of the
3 class:`cvxopt.base.matrix` class).
7 from cvxopt
import matrix
8 from cvxopt
.lapack
import gees
, gesdd
, syevr
13 def append_col(left
, right
):
15 Append two matrices side-by-side.
21 The two matrices to append to one another.
27 A new matrix consisting of ``right`` appended to the right
33 >>> A = matrix([1,2,3,4], (2,2))
34 >>> B = matrix([5,6,7,8,9,10], (2,3))
43 >>> print(append_col(A,B))
49 return matrix([left
.trans(), right
.trans()]).trans()
52 def append_row(top
, bottom
):
54 Append two matrices top-to-bottom.
60 The two matrices to append to one another.
66 A new matrix consisting of ``bottom`` appended below ``top``.
71 >>> A = matrix([1,2,3,4], (2,2))
72 >>> B = matrix([5,6,7,8,9,10], (3,2))
82 >>> print(append_row(A,B))
91 return matrix([top
, bottom
])
94 def eigenvalues(symmat
):
96 Return the eigenvalues of the given symmetric real matrix.
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.
106 The real symmetric matrix whose eigenvalues you want.
112 A list of the eigenvalues (in no particular order) of ``symmat``.
118 If the input matrix is not symmetric.
123 >>> A = matrix([[2,1],[1,2]], tc='d')
127 If the input matrix is not symmetric, it may not have real
128 eigenvalues, and we don't know what to do::
130 >>> A = matrix([[1,2],[3,4]])
132 Traceback (most recent call last):
134 TypeError: input must be a symmetric real matrix
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')
141 domain_dim
= symmat
.size
[0]
142 eigs
= matrix(0, (domain_dim
, 1), tc
='d')
144 # Create a copy of ``symmat`` here because ``syevr`` clobbers it.
145 dummy
= matrix(symmat
, symmat
.size
)
150 def eigenvalues_re(anymat
):
152 Return the real parts of the eigenvalues of the given square matrix.
158 The square matrix whose eigenvalues you want.
164 A list of the real parts (in no particular order) of the
165 eigenvalues of ``anymat``.
171 If the input matrix is not square.
176 This is symmetric and has two real eigenvalues:
178 >>> A = matrix([[2,1],[1,2]], tc='d')
179 >>> sorted(eigenvalues_re(A))
182 But this rotation matrix has eigenvalues `i` and `-i`, both of whose
185 >>> A = matrix([[0,-1],[1,0]])
186 >>> eigenvalues_re(A)
189 If the input matrix is not square, it doesn't have eigenvalues::
191 >>> A = matrix([[1,2],[3,4],[5,6]])
192 >>> eigenvalues_re(A)
193 Traceback (most recent call last):
195 TypeError: input matrix must be square
198 if not anymat
.size
[0] == anymat
.size
[1]:
199 raise TypeError('input matrix must be square')
201 domain_dim
= anymat
.size
[0]
202 eigs
= matrix(0, (domain_dim
, 1), tc
='z')
204 # Create a copy of ``anymat`` here for two reasons:
206 # 1. ``gees`` clobbers its input.
207 # 2. We need to ensure that the type code of ``dummy`` is 'd' or 'z'.
209 dummy
= matrix(anymat
, anymat
.size
, tc
='d')
212 return [eig
.real
for eig
in eigs
]
215 def identity(domain_dim
, typecode
='i'):
217 Create an identity matrix of the given dimensions.
223 The dimension of the vector space on which the identity will act.
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.
233 A ``domain_dim``-by-``domain_dim`` dense integer identity matrix.
239 If you ask for the identity on zero or fewer dimensions.
244 >>> print(identity(3))
252 raise ValueError('domain dimension must be positive')
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
)
260 def inner_product(vec1
, vec2
):
262 Compute the Euclidean inner product of two vectors.
268 The two vectors whose inner product you want.
274 The inner product of ``vec1`` and ``vec2``.
280 If the lengths of ``vec1`` and ``vec2`` differ.
287 >>> inner_product(x,y)
290 >>> x = matrix([1,1,1])
291 >>> y = matrix([2,3,4], (1,3))
292 >>> inner_product(x,y)
297 >>> inner_product(x,y)
298 Traceback (most recent call last):
300 TypeError: the lengths of vec1 and vec2 must match
303 if not len(vec1
) == len(vec2
):
304 raise TypeError('the lengths of vec1 and vec2 must match')
306 return sum([x
*y
for (x
, y
) in zip(vec1
, vec2
)])
309 def norm(matrix_or_vector
):
311 Return the Frobenius norm of a matrix or vector.
313 When the input is a vector, its matrix-Frobenius norm is the same
314 thing as its vector-Euclidean norm.
319 matrix_or_vector : matrix
320 The matrix or vector whose norm you want.
326 The norm of ``matrix_or_vector``.
331 >>> v = matrix([1,1])
332 >>> print('{:.5f}'.format(norm(v)))
335 >>> A = matrix([1,1,1,1], (2,2))
340 return sqrt(inner_product(matrix_or_vector
, matrix_or_vector
))
345 Create a long vector in column-major order from ``mat``.
351 Any sort of real matrix that you want written as a long vector.
357 An ``len(mat)``-by-``1`` long column vector containign the
358 entries of ``mat`` in column major order.
363 >>> A = matrix([[1,2],[3,4]])
376 Note that if ``mat`` is a vector, this function is a no-op:
378 >>> v = matrix([1,2,3,4], (4,1))
393 return matrix(mat
, (len(mat
), 1))
396 def condition_number(mat
):
398 Return the condition number of the given matrix.
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.
408 The larger the condition number is, the worse the matrix is.
413 The matrix whose condition number you want.
419 The nonnegative condition number of ``mat``.
424 >>> condition_number(identity(1, typecode='d'))
426 >>> condition_number(identity(2, typecode='d'))
428 >>> condition_number(identity(3, typecode='d'))
431 >>> A = matrix([[2,1],[1,2]], tc='d')
432 >>> abs(condition_number(A) - 3.0) < options.ABS_TOL
435 >>> A = matrix([[2,1j],[-1j,2]], tc='z')
436 >>> abs(condition_number(A) - 3.0) < options.ABS_TOL
440 num_eigs
= min(mat
.size
)
441 eigs
= matrix(0, (num_eigs
, 1), tc
='d')
445 return eigs
[0]/eigs
[-1]