X-Git-Url: http://gitweb.michael.orlitzky.com/?a=blobdiff_plain;f=src%2Fdunshire%2Fmatrices.py;h=ef3c4c0ba6b8926c67332b242a19c070032319ee;hb=2bf0987eea59afe20d1391babcd9b8562f7ec5b0;hp=0d97c0cbc4e83dfc95231a2306b870fb7b93a5db;hpb=814799cbabf73a2551532afa7f9e7128d9ae7bba;p=dunshire.git diff --git a/src/dunshire/matrices.py b/src/dunshire/matrices.py index 0d97c0c..ef3c4c0 100644 --- a/src/dunshire/matrices.py +++ b/src/dunshire/matrices.py @@ -1,20 +1,45 @@ """ Utility functions for working with CVXOPT matrices (instances of the -``cvxopt.base.matrix`` class). +class:`cvxopt.base.matrix` class). """ from math import sqrt from cvxopt import matrix -from cvxopt.lapack import syev +from cvxopt.lapack import gees, syevr + +import options + def append_col(left, right): """ - Append the matrix ``right`` to the right side of the matrix ``left``. + Append two matrices side-by-side. + + Parameters + ---------- + + left, right : matrix + The two matrices to append to one another. - EXAMPLES: + Returns + ------- + + matrix + A new matrix consisting of ``right`` appended to the right + of ``left``. + + Examples + -------- >>> A = matrix([1,2,3,4], (2,2)) >>> B = matrix([5,6,7,8,9,10], (2,3)) + >>> print(A) + [ 1 3] + [ 2 4] + + >>> print(B) + [ 5 7 9] + [ 6 8 10] + >>> print(append_col(A,B)) [ 1 3 5 7 9] [ 2 4 6 8 10] @@ -23,14 +48,37 @@ def append_col(left, right): """ return matrix([left.trans(), right.trans()]).trans() + def append_row(top, bottom): """ - Append the matrix ``bottom`` to the bottom of the matrix ``top``. + Append two matrices top-to-bottom. + + Parameters + ---------- + + top, bottom : matrix + The two matrices to append to one another. + + Returns + ------- + + matrix + A new matrix consisting of ``bottom`` appended below ``top``. - EXAMPLES: + Examples + -------- >>> A = matrix([1,2,3,4], (2,2)) >>> B = matrix([5,6,7,8,9,10], (3,2)) + >>> print(A) + [ 1 3] + [ 2 4] + + >>> print(B) + [ 5 8] + [ 6 9] + [ 7 10] + >>> print(append_row(A,B)) [ 1 3] [ 2 4] @@ -43,29 +91,148 @@ def append_row(top, bottom): return matrix([top, bottom]) -def eigenvalues(real_matrix): +def eigenvalues(symmat): """ - Return the eigenvalues of the given ``real_matrix``. + Return the eigenvalues of the given symmetric real matrix. + + On the surface, this appears redundant to the :func:`eigenvalues_re` + function. However, if we know in advance that our input is + symmetric, a better algorithm can be used. + + Parameters + ---------- + + symmat : matrix + The real symmetric matrix whose eigenvalues you want. + + Returns + ------- + + list of float + A list of the eigenvalues (in no particular order) of ``symmat``. + + Raises + ------ + + TypeError + If the input matrix is not symmetric. - EXAMPLES: + Examples + -------- >>> A = matrix([[2,1],[1,2]], tc='d') >>> eigenvalues(A) [1.0, 3.0] + If the input matrix is not symmetric, it may not have real + eigenvalues, and we don't know what to do:: + + >>> A = matrix([[1,2],[3,4]]) + >>> eigenvalues(A) + Traceback (most recent call last): + ... + TypeError: input must be a symmetric real matrix + """ - domain_dim = real_matrix.size[0] # Assume ``real_matrix`` is square. + if not norm(symmat.trans() - symmat) < options.ABS_TOL: + # Ensure that ``symmat`` is symmetric (and thus square). + raise TypeError('input must be a symmetric real matrix') + + domain_dim = symmat.size[0] eigs = matrix(0, (domain_dim, 1), tc='d') - syev(real_matrix, eigs) + syevr(symmat, eigs) return list(eigs) +def eigenvalues_re(anymat): + """ + Return the real parts of the eigenvalues of the given square matrix. + + Parameters + ---------- + + anymat : matrix + The square matrix whose eigenvalues you want. + + Returns + ------- + + list of float + A list of the real parts (in no particular order) of the + eigenvalues of ``anymat``. + + Raises + ------ + + TypeError + If the input matrix is not square. + + Examples + -------- + + This is symmetric and has two real eigenvalues: + + >>> A = matrix([[2,1],[1,2]], tc='d') + >>> sorted(eigenvalues_re(A)) + [1.0, 3.0] + + But this rotation matrix has eigenvalues `i` and `-i`, both of whose + real parts are zero: + + >>> A = matrix([[0,-1],[1,0]]) + >>> eigenvalues_re(A) + [0.0, 0.0] + + If the input matrix is not square, it doesn't have eigenvalues:: + + >>> A = matrix([[1,2],[3,4],[5,6]]) + >>> eigenvalues_re(A) + Traceback (most recent call last): + ... + TypeError: input matrix must be square + + """ + if not anymat.size[0] == anymat.size[1]: + raise TypeError('input matrix must be square') + + domain_dim = anymat.size[0] + eigs = matrix(0, (domain_dim, 1), tc='z') + + # Create a copy of ``anymat`` here for two reasons: + # + # 1. ``gees`` clobbers its input. + # 2. We need to ensure that the type code of ``dummy`` is 'd' or 'z'. + # + dummy = matrix(anymat, anymat.size, tc='d') + + gees(dummy, eigs) + return [eig.real for eig in eigs] + + def identity(domain_dim): """ - Return a ``domain_dim``-by-``domain_dim`` dense integer identity - matrix. + Create an identity matrix of the given dimensions. + + Parameters + ---------- + + domain_dim : int + The dimension of the vector space on which the identity will act. + + Returns + ------- - EXAMPLES: + matrix + A ``domain_dim``-by-``domain_dim`` dense integer identity matrix. + + Raises + ------ + + ValueError + If you ask for the identity on zero or fewer dimensions. + + Examples + -------- >>> print(identity(3)) [ 1 0 0] @@ -83,13 +250,76 @@ def identity(domain_dim): return matrix(entries, (domain_dim, domain_dim)) +def inner_product(vec1, vec2): + """ + Compute the Euclidean inner product of two vectors. + + Parameters + ---------- + + vec1, vec2 : matrix + The two vectors whose inner product you want. + + Returns + ------- + + float + The inner product of ``vec1`` and ``vec2``. + + Raises + ------ + + TypeError + If the lengths of ``vec1`` and ``vec2`` differ. + + Examples + -------- + + >>> x = [1,2,3] + >>> y = [3,4,1] + >>> inner_product(x,y) + 14 + + >>> x = matrix([1,1,1]) + >>> y = matrix([2,3,4], (1,3)) + >>> inner_product(x,y) + 9 + + >>> x = [1,2,3] + >>> y = [1,1] + >>> inner_product(x,y) + Traceback (most recent call last): + ... + TypeError: the lengths of vec1 and vec2 must match + + """ + if not len(vec1) == len(vec2): + raise TypeError('the lengths of vec1 and vec2 must match') + + return sum([x*y for (x, y) in zip(vec1, vec2)]) + + def norm(matrix_or_vector): """ - Return the Frobenius norm of ``matrix_or_vector``, which is the same - thing as its Euclidean norm when it's a vector (when one of its - dimensions is unity). + Return the Frobenius norm of a matrix or vector. + + When the input is a vector, its matrix-Frobenius norm is the same + thing as its vector-Euclidean norm. + + Parameters + ---------- + + matrix_or_vector : matrix + The matrix or vector whose norm you want. + + Returns + ------- + + float + The norm of ``matrix_or_vector``. - EXAMPLES: + Examples + -------- >>> v = matrix([1,1]) >>> print('{:.5f}'.format(norm(v))) @@ -100,4 +330,57 @@ def norm(matrix_or_vector): 2.0 """ - return sqrt(sum([x**2 for x in matrix_or_vector])) + return sqrt(inner_product(matrix_or_vector, matrix_or_vector)) + + +def vec(mat): + """ + Create a long vector in column-major order from ``mat``. + + Parameters + ---------- + + mat : matrix + Any sort of real matrix that you want written as a long vector. + + Returns + ------- + + matrix + An ``len(mat)``-by-``1`` long column vector containign the + entries of ``mat`` in column major order. + + Examples + -------- + + >>> A = matrix([[1,2],[3,4]]) + >>> print(A) + [ 1 3] + [ 2 4] + + + >>> print(vec(A)) + [ 1] + [ 2] + [ 3] + [ 4] + + + Note that if ``mat`` is a vector, this function is a no-op: + + >>> v = matrix([1,2,3,4], (4,1)) + >>> print(v) + [ 1] + [ 2] + [ 3] + [ 4] + + >>> print(vec(v)) + [ 1] + [ 2] + [ 3] + [ 4] + + + """ + return matrix(mat, (len(mat), 1))