1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
|
# Routines for creating and manipulating subspaces of vector spaces
# FIXME: Currently pretty broken, since instead of
# subspace objects, I have to work with matrices such that their
# column span is the subspace. (which i make sure are linearly indep.)
# FIXME: Also, I don't deal well with empty spaces/sets
# FIXME: Also, I do no parameter checking
## Vector Subspace creation
# Column Space of a matrix
#Column reduce and kill zero vectors
SetHelp ("ColumnSpace", "linear_algebra", "Get a basis matrix for the columnspace of a matrix")
function ColumnSpace(M) =
(
if IsNull (M) then return null
else if not IsMatrix (M) then
(error("ColumnSpace: argument not a matrix");bailout);
StripZeroColumns (cref (M))
)
# Row Space of a matrix
SetHelp ("RowSpace", "linear_algebra", "Get a basis matrix for the rowspace of a matrix")
function RowSpace(M) = (
if IsNull (M) then return null
else if not IsMatrix (M) then
(error("RowSpace: argument not a matrix");bailout);
ColumnSpace (M.')
)
# Rank of a matrix
SetHelp ("Rank", "linear_algebra", "Get the rank of a matrix")
function Rank(M) =
(
if IsNull (M) then return 0
else if not IsMatrix (M) then
(error("Rank: argument not a matrix");bailout);
columns (M) - CountZeroColumns(cref(M))
)
SetHelpAlias ("Rank", "rank");
rank = Rank
# Nullity of a matrix
SetHelp ("Nullity", "linear_algebra", "Get the nullity of a matrix")
function Nullity (M) =
(
if IsNull (M) then return 0
else if not IsMatrix (M) then
(error("Nullity: argument not a matrix");bailout);
CountZeroColumns(cref(M))
)
SetHelpAlias ("Nullity", "nullity");
nullity = Nullity
# Image of a linear transform (=column space of the corresponding matrix)
SetHelp ("Image", "linear_algebra", "Get the image (columnspace) of a linear transform")
Image = ColumnSpace
# Null space/kernel of a linear transform
# NullSpace is now built-in for speed
SetHelp ("Kernel", "linear_algebra", "Get the kernel (nullspace) of a linear transform")
function Kernel(M) = NullSpace(M)
## Vector Subspace operations
# Given W1, W2 subspaces of V, returns W1 + W2
# = {w | w=w1+w2, w1 in W1, w2 in W2}
SetHelp ("VectorSubspaceSum", "linear_algebra", "The sum of the vector spaces M and N, that is {w | w=m+n, m in M, n in N}")
function VectorSubspaceSum(M,N)=ColumnSpace([M,N])
# Given vector spaces V1, V2, return V1 (+) V2 (the direct sum)
# (given by the column space of the direct sum of matrices)
SetHelp ("VectorSpaceDirectSum", "linear_algebra", "The direct sum of the vector spaces M and N")
function VectorSpaceDirectSum(M,N)=ColumnSpace(DirectSum(M,N))
# Given vector spaces V1, V2, return V1 (x) V2 (the tensor product)
# (given by the column space of the tensor product of matrices)
#FIXME: check that this is right
#FIXME: need tensorproduct
# function VectorSpaceTensorProduct(M,N)=ColumnSpace(TensorProduct(M,N))
# Given W1, W2, return W1 intersect W2
# Given any inner product, this is given by
# W1 cap W2 = (W1^perp+W2^perp)^perp
SetHelp ("VectorSubspaceIntersection", "linear_algebra", "Intersection of the subspaces given by M and N")
function VectorSubspaceIntersection(M,N)= OrthogonalComplement(VectorSubspaceSum(OrthogonalComplement(M),OrthogonalComplement(N)))
# Given a vector subspace of an Inner Product Space (or maybe just
# a space with a (symmetric?) bilinear form (check lang),
# return the orthogonal complement.
# FIXME: currently assumes that the inner product is the hermetian
# inner product on
SetHelp ("OrthogonalComplement", "linear_algebra", "Get the orthogonal complement of the columnspace")
function OrthogonalComplement(M)=NullSpace(M')
## (Vector space + a vector) operations
# vector space membership
# a vector is in a subspace iff it is equal to its
# projection onto that space (in any inner product)
# (in particular, the dot product)
#FIXME: this is kinda a hack -- i should just ajoin it to a basis
#and see if there's a dependence (oh wait, that might not work for
#modules, since it might mean that a MULTIPLE of it is in the module
#....)
SetHelp ("IsInSubspace", "linear_algebra", "Test if a vector is in a subspace")
function IsInSubspace(v,W) = (v==Projection(v,W,.))
|
