Skip to content

jmichel7/GenLinearAlgebra.jl

Folders and files

NameName
Last commit message
Last commit date

Latest commit

 

History

12 Commits
 
 
 
 
 
 
 
 
 
 
 
 
 
 

Repository files navigation

Linear algebra on any field/ring

# GenLinearAlgebraModule.

GenLinearAlgebra: linear algebra over arbitrary fields and rings

The linear algebra package in Julia is not suitable for a general mathematics package: it assumes the field is the Real or Complex numbers and uses floating point to do approximate computations. For example the invertible rational matrix [1//(n+m) for n in 1:11, m in 1:11] has for LinearAlgebra a rank of 10 and a Float64 nullspace of dimension 1.

Here we are interested in functions which work over any field (or sometimes any ring) and assume exact computations.

For more information, look at the helpstrings of echelon!, rowspace, independent_rows, in_rowspace, intersect_rowspace, lnullspace, GenLinearAlgebra.nullspace, GenLinearAlgebra.rank, solutionmat, charpoly, comatrix, permanent, symmetric_power, exterior_power, ratio, diagconj_elt, transporter, bigcell_decomposition, traces_words_mats, all_ge_1.

source

# GenLinearAlgebra.echelon!Function.

echelon!(m::AbstractMatrix)

puts m in echelon form in-place and returns:

(m, indices of linearly independent rows of original m)

The echelon form transforms the rows of m into a particular basis of the rowspace. The first non-zero element of each line is 1, and such an element is also the only non-zero in its column. This function works in any field.

julia> m=[0 0 0 1; 1 1 0 1; 0 1 0 1; 1 0 0 0]//1
4×4 Matrix{Rational{Int64}}:
 0  0  0  1
 1  1  0  1
 0  1  0  1
 1  0  0  0

julia> echelon!(m)
(Rational{Int64}[1 0 0 0; 0 1 0 0; 0 0 0 1; 0 0 0 0], [2, 3, 1])

source

# GenLinearAlgebra.rowspaceFunction.

rowspace(m::AbstractMatrix)

returns a canonical form of the rowspace of m (the vector space generated by the rows of m), the echelon form. This is a particular basis of the rowspace of m: the first non-zero element of each line is 1, and such an element is also the only non-zero element in its column. Integer matrices are converted to Rational before applying this function.

julia> m=[1 2;2 4;5 6]
3×2 Matrix{Int64}:
 1  2
 2  4
 5  6

julia> rowspace(m)
2×2 view(::Matrix{Rational{Int64}}, 1:2, :) with eltype Rational{Int64}:
 1  0
 0  1

source

# GenLinearAlgebra.rankFunction.

GenLinearAlgebra.rank(m::AbstractMatrix)

computes exactly the rank of m. Not exported to not conflict with LinearAlgebra.

julia> hilbert(r)=[1//(n+m) for n in 1:r, m in 1:r];

julia> LinearAlgebra.rank(hilbert(11))
10

julia> GenLinearAlgebra.rank(hilbert(11)) # correct value
11

source

# GenLinearAlgebra.independent_rowsFunction.

independent_rows(m::AbstractMatrix)

returns a vector of the indices of a set of independent rows of m (the lexicographically first such set).

julia> m=[1 2;2 4;5 6]
3×2 Matrix{Int64}:
 1  2
 2  4
 5  6

julia> independent_rows(m)
2-element view(::Vector{Int64}, 1:2) with eltype Int64:
 1
 3

source

# GenLinearAlgebra.in_rowspaceFunction.

in_rowspace(v::AbstractVector,m::AbstractMatrix)

whether v is in the rowspace m. The matrix m should be an echelonized matrix like what is returned by function rowspace.

source

# GenLinearAlgebra.intersect_rowspaceFunction.

intersect_rowspace(m::AbstractMatrix,n::AbstractMatrix)

The intersection of the rowspaces of m and n.

source

# GenLinearAlgebra.lnullspaceFunction.

lnullspace(m::AbstractMatrix)

The left nullspace of m, that is the rowspace of the vectors v such that iszero(permutedims(v)*m).

source

# GenLinearAlgebra.nullspaceFunction.

GenLinearAlgebra.nullspace(m::AbstractMatrix)

computes the right nullspace of m in a type-preserving way, that is that is the colspace of the vectors v such that iszero(m*v). Not exported to avoid conflict with LinearAlgebra

source

# GenLinearAlgebra.solutionmatFunction.

solutionmat(mat,v)

returns one solution x of the equation transpose(mat)*x=v (that is, expresses v as a linear combination of the rows of mat), or nothing if no such solution exists. Similar to transpose(mat)\v when mat is invertible.

julia> m=[2 -4 1;0 0 -4;1 -2 -1]
3×3 Matrix{Int64}:
 2  -4   1
 0   0  -4
 1  -2  -1

julia> x=solutionmat(m,[10,-20, -10])
3-element Vector{Rational{Int64}}:
   5
 15//4
   0

julia> m'*x
3-element Vector{Rational{Int64}}:
  10
 -20
 -10

julia> solutionmat(m,[10, 20, -10])

source

solutionmat(m,n::AbstractMatrix)

return a matrix x such that x*m==n. This is interesting when m is not invertible.

source

# GenLinearAlgebra.charpolyFunction.

charpoly(M::Matrix)

The characteristic polynomial of M (as a Vector of coefficients). This function works over any ring.

source

# GenLinearAlgebra.comatrixFunction.

comatrix(M::Matrix)

is defined by comatrix(M)*M=det(M)*one(M). This function works over any ring.

source

# GenLinearAlgebra.permanentFunction.

permanent(m)

returns the permanent of the square matrix m, which is defined by $\sum_{p\in\frak S_n}\prod_{i=1}^n m[i,p(i)]$.

Note the similarity of the definition of the permanent to the definition of the determinant. In fact the only difference is the missing sign of the permutation. However the permanent is quite unlike the determinant, for example it is not multilinear or alternating. It has however important combinatorical properties.

julia> permanent([0 1 1 1;1 0 1 1;1 1 0 1;1 1 1 0]) # inefficient way to compute the number of derangements of 1:4
9

julia> permanent([1 1 0 1 0 0 0; 0 1 1 0 1 0 0;0 0 1 1 0 1 0; 0 0 0 1 1 0 1;1 0 0 0 1 1 0;0 1 0 0 0 1 1;1 0 1 0 0 0 1]) # 24 permutations fit the projective plane of order 2 
24 

source

# GenLinearAlgebra.symmetric_powerFunction.

symmetric_power(m,n)

returns the n-th symmetric power of the square matrix m, in the basis naturally indexed by the multisets of 1:n, where n=size(m,1).

julia> m=[1 2;3 4]
2×2 Matrix{Int64}:
 1  2
 3  4

julia> Int.(symmetric_power(m,2))
3×3 Matrix{Int64}:
 1   2   4
 6  10  16
 9  12  16

source

# GenLinearAlgebra.exterior_powerFunction.

exterior_power(mat,n)

mat should be a square matrix. The function returns the n-th exterior power of mat, in the basis naturally indexed bycombinations(1:r,n) wherer=size(mat,1)

julia> M=[1 2 3 4;2 3 4 1;3 4 1 2;4 1 2 3]
4×4 Matrix{Int64}:
 1  2  3  4
 2  3  4  1
 3  4  1  2
 4  1  2  3

julia> exterior_power(M,2)
6×6 Matrix{Int64}:
  -1   -2   -7   -1  -10  -13
  -2   -8  -10  -10  -12    2
  -7  -10  -13    1    2    1
  -1  -10    1  -13    2    7
 -10  -12    2    2    8   10
 -13    2    1    7   10   -1

source

# GenLinearAlgebra.ratioFunction.

ratio(v::abstractArray,w::abstractArray)

ratio r such that v=r.*w, nothing if v is not a multiple of w.

source

# GenLinearAlgebra.diagconj_eltFunction.

diagconj_elt(M,N)

M and N must be square matrices of the same size. This function returns a list d such that N==inv(Diagonal(d))*M*Diagonal(d) if such a list exists, and nothing otherwise.

julia> M=[1 2;2 1];N=[1 4;1 1]
2×2 Matrix{Int64}:
 1  4
 1  1

julia> diagconj_elt(M,N)
2-element Vector{Rational{Int64}}:
 1
 2

source

# GenLinearAlgebra.transporterFunction.

transporter(l1, l2 )

l1 and l2 should be vectors of the same length of square matrices all of the same size. The result is a basis of the vector space of matrices A such that for any i we have A*l1[i]=l2[i]*A –- the basis is returned as a vector of matrices, empty if the vector space is 0. This is useful to find whether two representations are isomorphic.

source

# GenLinearAlgebra.bigcell_decompositionFunction.

bigcell_decomposition(M [, b])

M should be a square matrix, and b specifies a block structure for a matrix of same size as M (it is a Vector of Vectors whose concatenation is 1:size(M,1)). If b is not given, the trivial block structure [[i] for i in axes(M,1)] is assumed.

The function decomposes M as a product P₁ L P where P is upper block-unitriangular (with identity diagonal blocks), P₁ is lower block-unitriangular and L is block-diagonal for the block structure b. If M is symmetric then P₁ is the transposed of P and the result is the pair [P,L]; else the result is the triple [P₁,L,P]. The only condition for this decomposition of M to be possible is that the principal minors according to the block structure be invertible. This routine is used in the Lusztig-Shoji algorithm for computing the Green functions and the example below is extracted from the computation of the Green functions for G₂.

julia> using LaurentPolynomials

julia> @Pol q
Pol{Int64}: q

julia> M=[q^6 q^0 q^3 q^3 q^5+q q^4+q^2; q^0 q^6 q^3 q^3 q^5+q q^4+q^2; q^3 q^3 q^6 q^0 q^4+q^2 q^5+q; q^3 q^3 q^0 q^6 q^4+q^2 q^5+q; q^5+q q^5+q q^4+q^2 q^4+q^2 q^6+q^4+q^2+1 q^5+2*q^3+q; q^4+q^2 q^4+q^2 q^5+q q^5+q q^5+2*q^3+q q^6+q^4+q^2+1]
6×6 Matrix{Pol{Int64}}:
 q⁶     1      q³     q³     q⁵+q        q⁴+q²
 1      q⁶     q³     q³     q⁵+q        q⁴+q²
 q³     q³     q⁶     1      q⁴+q²       q⁵+q
 q³     q³     1      q⁶     q⁴+q²       q⁵+q
 q⁵+q   q⁵+q   q⁴+q²  q⁴+q²  q⁶+q⁴+q²+1  q⁵+2q³+q
 q⁴+q²  q⁴+q²  q⁵+q   q⁵+q   q⁵+2q³+q    q⁶+q⁴+q²+1

julia> bb=[[2],[4],[6],[3,5],[1]];

julia> (P,L)=bigcell_decomposition(M,bb);

julia> P
6×6 Matrix{Pol{Int64}}:
 1    0  0    0    0        0
 q⁻⁶  1  q⁻³  q⁻³  q⁻¹+q⁻⁵  q⁻²+q⁻⁴
 0    0  1    0    0        0
 q⁻³  0  0    1    q⁻²      q⁻¹
 q⁻¹  0  0    0    1        0
 q⁻²  0  q⁻¹  0    q⁻¹      1

julia> L
6×6 Matrix{Pol{Int64}}:
 q⁶-q⁴-1+q⁻²  0   0            0     0            0
 0            q⁶  0            0     0            0
 0            0   q⁶-q⁴-1+q⁻²  0     0            0
 0            0   0            q⁶-1  0            0
 0            0   0            0     q⁶-q⁴-1+q⁻²  0
 0            0   0            0     0            q⁶-1

julia> M==transpose(P)*L*P
true

source

# GenLinearAlgebra.traces_words_matsFunction.

traces_words_mats(mats,words)

given a list mats of matrices and a list words of words returns the list of traces of the corresponding products of the matrices. Each subword is evaluated only once.

julia> R=[[-1 -1 0 0; 0 1 0 0; 0 0 1 0; 0 0 0 1], [1 0 0 0; -1 -1 -1 0; 0 0 1 0; 0 0 0 1], [1 0 0 0; 0 1 0 0; 0 -2 -1 -1; 0 0 0 1], [1 0 0 0; 0 1 0 0; 0 0 1 0; 0 0 -1 -1]]; # 17th representation of F4

julia> words=[[], [1,2,3,4,1,2,3,4,1,2,3,4,1,2,3,4,1,2,3,4], [2,3,2,3], [2,1], [1,2,3,4,2,3,2,3,4,3], [1,2,3,4,1,2,3,4,1,2,3,4], [4,3], [1,2,1,3,2,3,1,2,3,4], [1,2,3,4,1,2,3,4,1,2,3,4,1,2,3,4], [1,2,3,4,1,2,3,4], [1,2,3,4], [1], [2,3,2,3,4,3,2,3,4], [1,4,3], [4,3,2], [2,3,2,1,3], [3], [1,2,1,3,2,1,3,2,3], [2,1,4], [3,2,1], [2,4,3,2,3], [1,3], [3,2], [1,2,3,4,1,2,3,4,1,2,3,4,2,3], [1,2,3,4,2,3]];

julia> [traces_words_mats(R,words)] # 17th character of F4
1-element Vector{Vector{Int64}}:
 [4, 0, 0, 1, -1, 0, 1, -1, -2, 2  …  0, 2, -2, -1, 1, 0, 0, 2, -2, 0]

source

# GenLinearAlgebra.all_ge_1Function.

all_ge_1(M::Matrix;approx=x->x)

This function determines if a list of real linear forms represented by the rows of a matrix M can be made simultaneously ≥1. The result is a vector v such that all(≥(1),M*v), or nothing if there is no such vector.

approx(x) is a function giving an approximate real value for x. For instance, setting approx=Float64 enables to use this function with real Cyc numbers.

julia> all_ge_1([1 1 -1;1 -1 1;-1 1 1])
3-element Vector{Float64}:
 1.0
 1.0
 1.0

source

About

Linear algebra on an arbitrary field or ring

Resources

License

Stars

Watchers

Forks

Releases

No releases published

Packages

No packages published

Languages