CoCalc -- 2018-01-29.ipynb

| Download

This repository contains the course materials from Math 157: Intro to Mathematical Software.

Creative Commons BY-SA 4.0 license.

Project: Support and Testing

Path: math157 / lectures / 2018-01-29 / 2018-01-29.ipynb

Views: ³⁰³⁷
License: OTHER

Kernel: SageMath 8.1

Math 157: Intro to Mathematical Software

UC San Diego, winter 2018

January 29, 2018: Linear algebra (part 1 of 2)

Administrivia:

Homework 1 has been returned. Feedback has been added to your project. Grades are posted on TritonEd.
Homework 2 solutions have been distributed. We are expecting to have grades and feedback by Friday at noon.
Homework 3 is now posted. It is due Friday, February 2 at 8pm.
In addition, 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.

Added in class:

Some of you might find this documentary screening on Monday, February 5 of interest.
We are still working some kinks out of the script that computes participation grades. Let us know if you find any issues.

Vectors and matrices in Sage

A handy "cheat sheet" for matrices in Sage: the Sage Quick Reference for Linear Algebra. You might also find useful:

References for today's material:
- Linear algebra reference manual: http://doc.sagemath.org/html/en/reference/index.html#linear-algebra
- 1.5, 4.4, 4.16, of Sage For Undergraduates: http://www.gregorybard.com/Sage.html

In [90]:

v = vector((2, -1, 4))
w = vector([1, 0, -3])
print(v,w)

((2, -1, 4), (1, 0, -3))

In [92]:

v + 3*w # Addition, scalar multiplication

(5, -1, -5)

In [93]:

v.dot_product(w)

-10

In [94]:

v.cross_product(w)

(3, 10, 1)

In [95]:

v*w # Synonym for dot product

-10

In [97]:

v.norm(p=1)

7

In [0]:

v.norm?

In [98]:

M = Matrix([[1, -1, 0], [-2, 2, -1], [0, 1, 3]])
M.det() # Shorthand for M.determinant()

1

In [99]:

Matrix((v,w))

[ 2 -1  4]
[ 1  0 -3]

In [0]:

Matrix(2,2,[1,2,3,4])

In [101]:

M + identity_matrix(3)

[ 2 -1  0]
[-2  3 -1]
[ 0  1  4]

In [102]:

M * M # Matrix multiplication

[ 3 -3  1]
[-6  5 -5]
[-2  5  8]

In [103]:

M[1] # Extract a row

(-2, 2, -1)

In [104]:

M.transpose()[1] # Extract a column

(-1, 2, 1)

In [105]:

M.column(1)

(-1, 2, 1)

In [106]:

M[:2] # Extract a row slice, as a new matrix

[ 1 -1  0]
[-2  2 -1]

In [107]:

M.charpoly() # Short for M.characteristic_polynomial()

x^3 - 6*x^2 + 10*x - 1

Try this now: solve the system of linear equations $3x_1 + 2x_2 - x_3 = 1 \\ 2x_1 + x_2 - x_3 = -1 \\ x_1 - x_2 + x_3 = 2.$ There are several ways to do this! How many can you find?

In [114]:

# Try this here...
A = Matrix([[3,2,-1],[2,1,-1],[1,-1,1]])
v = vector([1,-1,2])
A.solve_right(v)

(1/3, 5/3, 10/3)

In [110]:

w = vector((1/3,5/3,10/3))
A*w

(1, -1, 2)

In [112]:

A.inverse()*v

(1/3, 5/3, 10/3)

In [118]:

l = A.columns()
l.append(v)
Matrix(l).transpose()

[ 3  2 -1  1]
[ 2  1 -1 -1]
[ 1 -1  1  2]

In [119]:

A.augment(v)

[ 3  2 -1  1]
[ 2  1 -1 -1]
[ 1 -1  1  2]

In [113]:

B = Matrix([[3,2,-1,1],[2,1,-1,-1],[1,-1,1,2]])
B.rref()

[   1    0    0  1/3]
[   0    1    0  5/3]
[   0    0    1 10/3]

Matrix entries

Unlike a Python array, the entries of a matrix must all be of the same type. Moreover, this type must admit addition, subtraction, and multiplication (that is, it must be an element of a ring). For some operations (e.g., row reduction), one also needs division by nonzero elements (i.e., a field).

Examples:

Integers (parent ring Integers() == IntegerRing() == ZZ)
Rational numbers (parent ring Rationals() == QQ)
Real numbers (parent ring RealField() == RR)
Complex numbers (parent ring ComplexField() == CC)
Can also change precision: RealField(200) for 200 bits of precision
Interval arithmetic: RealIntervalField(200), ComplexIntervalField(200)
Integers modulo n (parent ring Integers(n) == IntegerModRing(n))
For integers modulo a prime p, another synonym for the parent ring is GF(p)
Symbolic expressions (parent ring SymbolicRing() == SR)

To pass from one ring to another, use the change_ring function.

In [120]:

M = Matrix([[1, -1, 0], [-2, 2, -1], [0, 1, 3]])
N = M.change_ring(GF(2)) # Entries are now integers mod 2
N.inverse()

[1 1 1]
[0 1 1]
[0 1 0]

In [121]:

M = Matrix([[sqrt(2), pi], [e, I]])
M.inverse()

[-pi*e/(sqrt(2)*pi*e - 2*I) + 1/2*sqrt(2)          sqrt(2)*pi/(sqrt(2)*pi*e - 2*I)]
[          sqrt(2)*e/(sqrt(2)*pi*e - 2*I)                  -2/(sqrt(2)*pi*e - 2*I)]

In [122]:

N = M.change_ring(CC) # Can't use RR here because of the I in the matrix
N.inverse()

[ 0.0188745367983670 - 0.113973965547083*I  0.358059772863211 + 0.0592961061456607*I]
[ 0.309813359464052 + 0.0513063103995832*I -0.161183527834128 - 0.0266926259237610*I]

Constructing a matrix from a formula

The matrix constructor accepts inputs of several different types. One particularly useful syntax is to specify the dimensions (rows, then columns) and then a function that takes a row and column index (which as usual start from 0) and returns the corresponding matrix entry. This function can either be an explicitly declared function or a lambda function.

In [123]:

def f(i,j):
    return binomial(i,j)

M = Matrix(10, 10, f); M

[  1   0   0   0   0   0   0   0   0   0]
[  1   1   0   0   0   0   0   0   0   0]
[  1   2   1   0   0   0   0   0   0   0]
[  1   3   3   1   0   0   0   0   0   0]
[  1   4   6   4   1   0   0   0   0   0]
[  1   5  10  10   5   1   0   0   0   0]
[  1   6  15  20  15   6   1   0   0   0]
[  1   7  21  35  35  21   7   1   0   0]
[  1   8  28  56  70  56  28   8   1   0]
[  1   9  36  84 126 126  84  36   9   1]

In [124]:

M = Matrix(10, 10, lambda i,j: binomial(i, j))
M

[  1   0   0   0   0   0   0   0   0   0]
[  1   1   0   0   0   0   0   0   0   0]
[  1   2   1   0   0   0   0   0   0   0]
[  1   3   3   1   0   0   0   0   0   0]
[  1   4   6   4   1   0   0   0   0   0]
[  1   5  10  10   5   1   0   0   0   0]
[  1   6  15  20  15   6   1   0   0   0]
[  1   7  21  35  35  21   7   1   0   0]
[  1   8  28  56  70  56  28   8   1   0]
[  1   9  36  84 126 126  84  36   9   1]

In [125]:

M.inverse()

[   1    0    0    0    0    0    0    0    0    0]
[  -1    1    0    0    0    0    0    0    0    0]
[   1   -2    1    0    0    0    0    0    0    0]
[  -1    3   -3    1    0    0    0    0    0    0]
[   1   -4    6   -4    1    0    0    0    0    0]
[  -1    5  -10   10   -5    1    0    0    0    0]
[   1   -6   15  -20   15   -6    1    0    0    0]
[  -1    7  -21   35  -35   21   -7    1    0    0]
[   1   -8   28  -56   70  -56   28   -8    1    0]
[  -1    9  -36   84 -126  126  -84   36   -9    1]

Try this on your own: formulate a general conjecture based on this example, then prove it. Better yet, try variations of this construction to see if you can come with any other conjectures like it!

Numpy

The NumPy ("Numerical Python") project provides Python with functionality comparable to Fortran. From the home page:

NumPy is the fundamental package for scientific computing with Python. It contains among other things:
a powerful N-dimensional array object
sophisticated (broadcasting) functions
tools for integrating C/C++ and Fortran code
useful linear algebra, Fourier transform, and random number capabilities

NumPy is a pure Python package, but requires installation. However, Sage provides NumPy by default.

In [126]:

import numpy as np # A handy shorthand for importing modules
a = np.array([[2,3,4],[5,6,7]]); a

array([[2, 3, 4],
       [5, 6, 7]])

The NumPy type array plays the role of both vectors and matrices. In fact, arrays may have any number of dimensions.

In [127]:

a = np.ones((2,3,4)); a

array([[[ 1.,  1.,  1.,  1.],
        [ 1.,  1.,  1.,  1.],
        [ 1.,  1.,  1.,  1.]],

       [[ 1.,  1.,  1.,  1.],
        [ 1.,  1.,  1.,  1.],
        [ 1.,  1.,  1.,  1.]]])

Converting between Sage matrices and Numpy matrices is easy! This makes it possible to apply numpy functions to Sage matrices.

In [128]:

M = Matrix([[1, -1, 0], [-2, 2, -1], [0, 1, 3]])
a = M.numpy(); a

array([[ 1, -1,  0],
       [-2,  2, -1],
       [ 0,  1,  3]])

Here we use the NumPy cumulative sum function.

In [131]:

b = np.cumsum(a, 0); b # Cumulative sum function

array([[ 1, -1,  0],
       [-1,  1, -1],
       [-1,  2,  2]])

In [132]:

Matrix(b)

[ 1 -1  0]
[-1  1 -1]
[-1  2  2]

Some other examples of NumPy functions you might want to access from Sage:

least squares approximations (see homework);
matrix decompositions (e.g., singular value decomposition);
discrete Fourier transform.

NumPy and SciPy

A closely related project to NumPy is SciPy ("Scientific Python"). More precisely, I am talking about the SciPy library, which together with NumPy form two components of the SciPy ecosystem. (Other components of the ecosystem include IPython, the basis for the Jupyter notebook; Matplotlib and Sympy, used by Sage for plotting and symbolics; and pandas, more on which later.)

Some of the strengths of SciPy (the library):

low-rank matrix approximations (a/k/a compressed sensing);
interpolation and curve-fitting;
numerical integration (including numerical solution of differential equations);
image manipulation;
signal processing.

We will probably not cover much of this in this course, though. See the SciPy lectures for a number of worked examples.

In [137]:

# Example of numerical solution of an ODE from SciPy lectures 1.5.7.2
# y' = -2y

def calc_derivative(ypos, time):
    return -2 * ypos

n = 400
from scipy.integrate import odeint
time_vec = np.linspace(0, 4, n)
y = odeint(calc_derivative, y0=1, t=time_vec)

In [134]:

time_vec

array([ 0.        ,  0.1025641 ,  0.20512821,  0.30769231,  0.41025641,
        0.51282051,  0.61538462,  0.71794872,  0.82051282,  0.92307692,
        1.02564103,  1.12820513,  1.23076923,  1.33333333,  1.43589744,
        1.53846154,  1.64102564,  1.74358974,  1.84615385,  1.94871795,
        2.05128205,  2.15384615,  2.25641026,  2.35897436,  2.46153846,
        2.56410256,  2.66666667,  2.76923077,  2.87179487,  2.97435897,
        3.07692308,  3.17948718,  3.28205128,  3.38461538,  3.48717949,
        3.58974359,  3.69230769,  3.79487179,  3.8974359 ,  4.        ])

In [135]:

y #Let's see what we got here.

array([[  1.00000000e+00],
       [  8.14542898e-01],
       [  6.63480126e-01],
       [  5.40433015e-01],
       [  4.40205868e-01],
       [  3.58566566e-01],
       [  2.92067826e-01],
       [  2.37901766e-01],
       [  1.93781189e-01],
       [  1.57843082e-01],
       [  1.28569955e-01],
       [  1.04725738e-01],
       [  8.53036015e-02],
       [  6.94834394e-02],
       [  5.65972376e-02],
       [  4.61008805e-02],
       [  3.75511484e-02],
       [  3.05870204e-02],
       [  2.49144403e-02],
       [  2.02938795e-02],
       [  1.65302356e-02],
       [  1.34645860e-02],
       [  1.09674830e-02],
       [  8.93348546e-03],
       [  7.27670714e-03],
       [  5.92719016e-03],
       [  4.82795066e-03],
       [  3.93257294e-03],
       [  3.20324937e-03],
       [  2.60918427e-03],
       [  2.12529267e-03],
       [  1.73114166e-03],
       [  1.41008843e-03],
       [  1.14857711e-03],
       [  9.35565007e-04],
       [  7.62057552e-04],
       [  6.20728409e-04],
       [  5.05609460e-04],
       [  4.11840447e-04],
       [  3.35461687e-04]])

In [138]:

# Collate the x and y axes and make a list plot.
point([(time_vec[i], y[i]) for i in range(n)])

In [0]: