X-Git-Url: http://gitweb.michael.orlitzky.com/?a=blobdiff_plain;f=src%2Fdunshire%2Fmatrices.py;h=38d4afe78efaad7665b71d936c8ba9e8309de581;hb=479cd91a21bb68b88021575741d812378971a7c9;hp=1e0f2a5b7970fd19b78a2d136dc0a98dee2547af;hpb=002b5370da24f083d2088c3482cf076615a13563;p=dunshire.git diff --git a/src/dunshire/matrices.py b/src/dunshire/matrices.py index 1e0f2a5..38d4afe 100644 --- a/src/dunshire/matrices.py +++ b/src/dunshire/matrices.py @@ -1,20 +1,36 @@ """ 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 +import options + def append_col(left, right): """ Append the matrix ``right`` to the right side of the matrix ``left``. - EXAMPLES: + Parameters + ---------- + left, right : matrix + The two matrices to append to one another. + + 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] @@ -27,10 +43,25 @@ def append_row(top, bottom): """ Append the matrix ``bottom`` to the bottom of the matrix ``top``. - EXAMPLES: + Parameters + ---------- + top, bottom : matrix + The two matrices to append to one another. + + 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 +74,66 @@ 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. + + Parameters + ---------- + + symmat : matrix + The real symmetric matrix whose eigenvalues you want. + - EXAMPLES: + Raises + ------ + TypeError + If the input matrix is not symmetric. + + 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) + syev(symmat, eigs) return list(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. - EXAMPLES: + Returns + ------- + matrix + A ``domain_dim``-by-``domain_dim`` dense integer identity matrix. + + Examples + -------- >>> print(identity(3)) [ 1 0 0] @@ -83,13 +151,69 @@ 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``. + + 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,27 +224,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(real_matrix): +def vec(mat): """ - Create a long vector in column-major order from ``real_matrix``. + Create a long vector in column-major order from ``mat``. - EXAMPLES: + Parameters + ---------- - >>> A = matrix([[1,2],[3,4]]) - >>> print(A) - [ 1 3] - [ 2 4] - + mat : matrix + Any sort of real matrix that you want written as a long vector. - >>> print(vec(A)) - [ 1] - [ 2] - [ 3] - [ 4] - + 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(real_matrix, (len(real_matrix), 1)) + return matrix(mat, (len(mat), 1))