2018-01-31.ipynb

Kernel: SageMath 8.1

January 31, 2018: Linear algebra (part 2 of 2)

Math 157: Intro to Mathematical Software

UC San Diego, winter 2018

Administrivia:

• Please fill out the week 4 questionnaire.

• We are continuing to monitor the course waitlist. If you are on it and still want to join the course, please continue to keep up with lectures and homework, and watch your email for further instructions.

• Here is a more precise link for the documentary screening next Monday that I mentioned last time.

Octave from Sage

So far we have seen Octave and Sage being accessed from a Jupyter notebook by switching kernels back and forth. However, this approach does not allow for any communication between the two kernels: all state is lost when one switches the kernel from Sage to Octave or vice versa.

It was asked whether there is a way to combine the functionality of the two kernels. The answer is yes, in one direction: Sage can call into Octave, but not vice versa. Some examples from the reference manual:

octave("besselj(1,2)")      # directly call a single Octave command

 0.576725

my_vector1 = octave('[1,5,7]') # create Python/Sage variables which point to Octave objects
my_vector2 = octave('[1;5;7]') 
my_vector1 * my_vector2 # Sage translates this operation into Octave for you.

 75

M33 = MatrixSpace(QQ,3,3)
A   = M33([1,2,3,4,5,6,7,8,0])
V3  = VectorSpace(QQ,3)
b   = V3([1,2,3])
octave.solve_linear_system(A,b)

[-0.333333, 0.666667, 0]

Vector spaces

So far we have been talking about linear algebra solely in terms of matrices. However, Sage can also handle the language of vector spaces, as long as we are talking about subspaces of a "standard" vector space over some field.

V = VectorSpace(QQ, 5)
V

Vector space of dimension 5 over Rational Field

v = V([2,1,1,-1,1])
print(v)
print(type(v))

(2, 1, 1, -1, 1)
<type 'sage.modules.vector_rational_dense.Vector_rational_dense'>

type(vector([2,1,1,-1,1]))

<type 'sage.modules.vector_integer_dense.Vector_integer_dense'>

w = V([1,1,-1,0,2])
W = V.span([v,w]); W

Vector space of degree 5 and dimension 2 over Rational Field
Basis matrix:
[ 1  0  2 -1 -1]
[ 0  1 -3  1  3]

W.ambient_vector_space() == V

True

w2 = V([1,1,-1,0,2])
w3 = V([2,-2,0,5,-1])
W2 = V.span([w2, w3])

W.intersection(W2)

Vector space of degree 5 and dimension 1 over Rational Field
Basis matrix:
[ 1  1 -1  0  2]

M = Matrix(QQ, [[2,3,1,4],[0,1,-1,2]])
M.right_kernel()

Vector space of degree 4 and dimension 2 over Rational Field
Basis matrix:
[   1    0 -2/3 -1/3]
[   0    1 -1/3 -2/3]

M = Matrix([[2,0,0],[0,3,1],[-1,1,2]])
M.eigenspaces_right()

[
(2, Vector space of degree 3 and dimension 1 over Rational Field
User basis matrix:
[ 1  1 -1]),
(1.381966011250106?, Vector space of degree 3 and dimension 1 over Algebraic Field
User basis matrix:
[                  0                   1 -1.618033988749895?]),
(3.618033988749895?, Vector space of degree 3 and dimension 1 over Algebraic Field
User basis matrix:
[                 0                  1 0.618033988749895?])
]

M = Matrix([[0,0],[1,0]])
M.eigenspaces_right()

[
(0, Vector space of degree 2 and dimension 1 over Rational Field
User basis matrix:
[0 1])
]

A = matrix(QQ, 3, [1..9])
show(A)
V = A.image()
W = A.kernel()
print "image: "
print V
print "kernel: "
print W

image: 
Vector space of degree 3 and dimension 2 over Rational Field
Basis matrix:
[ 1  0 -1]
[ 0  1  2]
kernel: 
Vector space of degree 3 and dimension 1 over Rational Field
Basis matrix:
[ 1 -2  1]

A.restrict(V) # V = image

[-6 -6]
[18 21]

A.restrict(W)

[0]

show(factor(A.charpoly()))

show(factor(A.restrict(V).charpoly()))

M = span(QQ, [vector([1,0,0])]) 
M
A.restrict(M)   # A doesn't leave the subspace M invariant... so this will fail!

---------------------------------------------------------------------------
ArithmeticError                           Traceback (most recent call last)
<ipython-input-42-eff670de8383> in <module>()
      1 M = span(QQ, [vector([Integer(1),Integer(0),Integer(0)])])
      2 M
----> 3 A.restrict(M)   # A doesn't leave the subspace M invariant... so this will fail!

/ext/sage/sage-8.1/src/sage/matrix/matrix2.pyx in sage.matrix.matrix2.Matrix.restrict (build/cythonized/sage/matrix/matrix2.c:41441)()
   5056                 C = [V.coordinate_vector(b*self) for b in V.basis()]
   5057             except ArithmeticError:
-> 5058                 raise ArithmeticError("subspace is not invariant under matrix")
   5059             return self.new_matrix(n, n, C, sparse=False)
   5060 
ArithmeticError: subspace is not invariant under matrix

Sparse vs. dense matrices

In computational linear algebra, it is important to distinguish between sparse matrices, in which all but a few entries are zero, and dense matrices, where most of the entries are nonzero. This is not a formal mathematical distinction among matrices, but rather a distinction between two different classes of algorithms for handling matrices.

v = vector({2:4, 95:4, 1000000:3/4})
parent(v)

Sparse vector space of dimension 1000001 over Rational Field

v[3]

0

M = Matrix({(3,2): 27, (2,98): 71})
parent(M)

Full MatrixSpace of 4 by 99 sparse matrices over Integer Ring

print(M.is_sparse(), M.is_dense())

(True, False)

Note that in the previous examples, Sage tried to guess the dimension of the matrix. If it is guessing wrong, you may need to fix this by specifying dimensions.

M = Matrix(QQ, 1000001, 1000001, {(3,2): 27, (2,98): 71})
parent(M)

Full MatrixSpace of 1000001 by 1000001 sparse matrices over Rational Field

M2 = M*M # Do not try this with dense matrices of this size!

w = M*v # Don't try to print this answer either.

sum(i for i in w)

Warning: Not all operations make sense to do in the sparse model. For example, passing from a matrix to its inverse does not preserve sparsity in general. Many operations on sparse matrices involve an implicit coercion to a dense matrix.

Operations on matrices

A = matrix(QQ, [[1,2],[3,4]])
A^3 - 2*A + 1

[ 36  50]
[ 75 111]

A+1

[2 2]
[3 5]

A = matrix(QQ, [[1,2],[3,4]])
f(x) = x^3 - 2*x + 1
f(A) # FAIL -- this should work, but doesn't

---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
<ipython-input-3-ce7578e55d8d> in <module>()
      1 A = matrix(QQ, [[Integer(1),Integer(2)],[Integer(3),Integer(4)]])
      2 __tmp__=var("x"); f = symbolic_expression(x**Integer(3) - Integer(2)*x + Integer(1)).function(x)
----> 3 f(A) # FAIL -- this should work, but doesn't

/ext/sage/sage-8.1/src/sage/symbolic/expression.pyx in sage.symbolic.expression.Expression.__call__ (build/cythonized/sage/symbolic/expression.cpp:32775)()
   5273             z^2 + x^y
   5274         """
-> 5275         return self._parent._call_element_(self, *args, **kwds)
   5276 
   5277     def variables(self):
/ext/sage/sage-8.1/local/lib/python2.7/site-packages/sage/symbolic/callable.pyc in _call_element_(self, _the_element, *args, **kwds)
    462         d = dict(zip([repr(_) for _ in self.arguments()], args))
    463         d.update(kwds)
--> 464         return SR(_the_element.substitute(**d))
    465 
    466     # __reduce__ gets replaced by the CallableSymbolicExpressionRingFactory
/ext/sage/sage-8.1/src/sage/symbolic/expression.pyx in sage.symbolic.expression.Expression.substitute (build/cythonized/sage/symbolic/expression.cpp:32085)()
   5171         for k, v in sdict.iteritems():
   5172             smap.insert(make_pair((<Expression>self.coerce_in(k))._gobj,
-> 5173                                   (<Expression>self.coerce_in(v))._gobj))
   5174         sig_on()
   5175         try:
/ext/sage/sage-8.1/src/sage/symbolic/expression.pyx in sage.symbolic.expression.Expression.coerce_in (build/cythonized/sage/symbolic/expression.cpp:23189)()
   3042             return <Expression?>z
   3043         except TypeError:
-> 3044             return self._parent._coerce_(z)
   3045 
   3046     cpdef _add_(left, right):
/ext/sage/sage-8.1/src/sage/structure/parent_old.pyx in sage.structure.parent_old.Parent._coerce_ (build/cythonized/sage/structure/parent_old.c:4962)()
    228     def _coerce_(self, x):            # Call this from Python (do not override!)
    229         if self._element_constructor is not None:
--> 230             return self.coerce(x)
    231         check_old_coerce(self)
    232         return self._coerce_c(x)
/ext/sage/sage-8.1/src/sage/structure/parent.pyx in sage.structure.parent.Parent.coerce (build/cythonized/sage/structure/parent.c:10941)()
   1178                 except Exception:
   1179                     _record_exception()
-> 1180             raise TypeError("no canonical coercion from %s to %s" % (parent(x), self))
   1181         else:
   1182             return (<map.Map>mor)._call_(x)
TypeError: no canonical coercion from Full MatrixSpace of 2 by 2 dense matrices over Rational Field to Callable function ring with argument x

R.<y> = PolynomialRing(QQ)
parent(y)

Univariate Polynomial Ring in y over Rational Field

A = matrix(QQ, [[1,2],[3,4]])
R.<y> = PolynomialRing(QQ)
f = y^3 - 2*y + 1
f(A)  # works

[ 36  50]
[ 75 111]

show(A.exp())   # compute sum A^k/k!
type(A.exp())

<type 'sage.matrix.matrix_symbolic_dense.Matrix_symbolic_dense'>

RR

Real Field with 53 bits of precision

RDF

Real Double Field

A.change_ring(RDF).exp()

[51.968956198705044  74.73656456700328]
[112.10484685050491 164.07380304920997]

print("A.inverse()")
print(A.inverse())
print("A^(-1)")
print(A^(-1))
print("~A")
print(~A)
print("A.transpose()")
print(A.transpose())
print("A.conjugate() -- boring since real")
print(A.conjugate())
print("A.antitranspose()")
print(A.antitranspose())
print("A.adjoint()")
print(A.adjoint())

A.inverse()
[  -2    1]
[ 3/2 -1/2]
A^(-1)
[  -2    1]
[ 3/2 -1/2]
~A
[  -2    1]
[ 3/2 -1/2]
A.transpose()
[1 3]
[2 4]
A.conjugate() -- boring since real
[1 2]
[3 4]
A.antitranspose()
[4 2]
[3 1]
A.adjoint()
[ 4 -2]
[-3  1]

More on NumPy

import numpy as np

v = np.arange(15) 
v

array([ 0,  1,  2,  3,  4,  5,  6,  7,  8,  9, 10, 11, 12, 13, 14])

a = v.reshape(3,5)
a

array([[ 0,  1,  2,  3,  4],
       [ 5,  6,  7,  8,  9],
       [10, 11, 12, 13, 14]])

print(a.shape, a.ndim, a.dtype)

((3, 5), 2, dtype('int64'))

Notice that I typed a.shape rather than a.shape() as you may have expected if you have gotten accustomed to Python syntax. The point is that shape is not a method of a NumPy array but rather an attribute, essentially a private variable belonging to the array.

c = np.array([ [1,2], [3,4] ], dtype=complex) # The entry type can be forced at creation
c

array([[ 1.+0.j,  2.+0.j],
       [ 3.+0.j,  4.+0.j]])

The following is an example of a vectorized operation, where we apply a function to an array, which applies that function componentwise. NumPy is set up to do this sort of thing very efficiently (a la Fortran).

v = np.linspace(0, 2, 5); v

array([ 0. ,  0.5,  1. ,  1.5,  2. ])

w = np.sin(v); w

array([ 0.        ,  0.47942554,  0.84147098,  0.99749499,  0.90929743])

Let's test out the speed using the Python function timeit.

v = np.arange(0,1,1000000)
timeit('w = np.sin(v)')

625 loops, best of 3: 1.69 µs per loop

timeit('w = np.array([np.sin(x) for x in v])')

625 loops, best of 3: 3.69 µs per loop

Observe the difference:

print("Numpy:")
a = np.array([[1,2],[3,4]])
print(np.exp(a))

Numpy:
[[  2.71828183   7.3890561 ]
 [ 20.08553692  54.59815003]]

print("Sage:")
b = matrix(RDF,[[1,2],[3,4]]) # RDF = RealDoubleField()
print(exp(b))

Sage:
[51.968956198705044  74.73656456700328]
[112.10484685050491 164.07380304920997]

print("Numpy:")
print(a * a)
print("Sage:")
print(b * b)

Numpy:
[[ 1  4]
 [ 9 16]]
Sage:
[ 7.0 10.0]
[15.0 22.0]

This is just an arbitrary choice, made differently by mathematicians and engineers, so you have to be aware of it.

np.exp(a)

array([[  2.71828183,   7.3890561 ],
       [ 20.08553692,  54.59815003]])

b.apply_map(exp)

[ 2.718281828459045   7.38905609893065]
[20.085536923187668 54.598150033144236]

b.apply_map(lambda x: x**2)

[ 1.0  4.0]
[ 9.0 16.0]

# and...
a.dot(a)

array([[ 7, 10],
       [15, 22]])

Some fun

matrix_plot(random_matrix(RDF,50)^3, cmap='summer', colorbar=True)

matrix_plot(matrix(50,50,lambda i,j: randint(0,1)))

v = vector([1,2,3])
w = vector([0,2,5])

var('s, t')
g = parametric_plot3d(s*v + t*w, (s, -1,1), (t,-1,1), opacity=.5)

v = vector([1,-2,3])
w = vector([0,3,-5])
h = parametric_plot3d(s*v + t*w, (s, -1,1), (t,-1,1), color='red')

#i = parametric_plot3d(...., (s, -2, 2), ...)
#g+h+i
show(g+h)