# Determinant and Jacobian

Today I started my linear algebra class by motivating the determinant using a little bit of history and applications. In particular I talked about determinant appearing in the change of variables in integration of multivariate functions in calculus:Then after the class I noticed that my picture of what happens when we linearize the image wasn’t very clear as I had drawn it. So I decided to take a look at an actual example. The following code takes a function $f$ that maps $\mathbb{R}^2$ to $\mathbb{R}^2$ and draws three regions:

• the unit cube in blue,
• the image of the unit cube under $f$ in red, and
• the linearization of the image by Jacobian of $f$ at the origin
def show_linearization(f,n=100):
# f is a function that maps R^2 to R^2 for example
# f(x,y) = ((x^2+2*(x+2)*(y+1))/4-1,y^2-(x+1)*y)
# The output is three areas:
# the unit cube
# the image of unit cube under f
# the linearization of the image using the Jacobian of f at (0,0)

var('x y')
A = jacobian(f,(x,y))(0,0)
p = f(0,0)
domxy = []
imgxy = []
jacxy = []
for i in range(n+1):
for j in range(n+1):
domxy.append((i/n,j/n))
imgxy.append(f(i/n,j/n))
jacxy.append(p+A*vector([i/n,j/n]))

P = points(domxy,color="blue",aspect_ratio=1, alpha = .3)
Q = points(imgxy,color="red",aspect_ratio=1, alpha = .3)
R = points(jacxy,color="green",aspect_ratio=1, alpha = .3)

(P+Q+R).show()

Here is a sample run:

f(x,y) = ((x^2+2*(x+2)*(y+1))/4-1,y^2-(x+1)*y)
show_linearization(f)

and the output is:

I chose the function in a way that $f(0,0)=(0,0)$ for simplicity. Here it is on sage cell server for you to play around with it: link.

What do you think? Any suggestions or comments?

# Linear approximation and Taylor polynomials

For a given function $f(x)$ around a point $a$, if you zoom in close enough on the graph of the function, you can’t tell the difference between the curve and the tangent line to the graph of the function at point $a$. This is the basic idea for linear approximations. To explore that statement here I wrote a little sage interactive code that will let you play around with it:

@interact
def _(
f = ('$f(x)=$', input_box(x^3 + 6*x^2 - 15*x + 7,width=20)),
a = ('around the point $x = a$', input_box(0,width=10)),
zoom = ('zoom', slider(.1,12,default=1))):

g(x) = f(a) + diff(f,x)(a)*(x-a)
P = plot(f,(x,a-10*2^(-zoom), a+10*2^(-zoom)),  color="blue", axes=True, thickness=4) # for horizontal tangent lines add 'aspect_ratio=1' to see something meaningful, otherwise sage will just change the ration so that the graph is not just a horizontal line.
Q = plot(g,(x,a-10*2^(-zoom), a+10*2^(-zoom)),  color="red", axes=True, thickness=2)
show(P + Q)


Click here to play with it in the sage cell server. Here is an output:

The linear approximations have errors though. One way to get better approximations is by taking into account the higher order derivatives of the function at the given point. The following code will calculate and draw the n-th degree Taylor polynomial of the function so that you can compare them:

@interact
def _(
f = ('$f(x)=$', input_box(sin(x),width=20)),
n = ('Degree of the Taylor polynomial, $n=$', slider([0..9],default=3)),
a = ('around the point $x = a$', input_box(0,width=10)),
xzoom = ('range of $x$', range_slider(-12,12,default=(-12,12))),
yzoom = ('range of $y$', range_slider(-20,20,default=(-4,4)))):

g = f.taylor(x, a, n)
P = plot(f,(x,xzoom[0], xzoom[1]), ymin = yzoom[0], ymax = yzoom[1], color="blue", axes=True, aspect_ratio=1, thickness=3)
Q = plot(g,(x,xzoom[0], xzoom[1]), ymin = yzoom[0], ymax = yzoom[1], color="red", axes=True, aspect_ratio=1)
R = text(r"$%s \approx %s$" %(f,g),(0,yzoom[1]+1),rgbcolor=(1,0,0))
show(P + Q + R)


Click here to see it in the sage cell server. Here is an output:

The only problem with that code is that it doesn’t render the latex formula correctly. That is, the two digit exponents are shown wrong. So, for now I’ve limited $n$ to 9. If you have any suggestions on how to fix it please see the discussion here.

# Some March Mathness

I went to several talks and meetings in the past few days, and I am going to talk about some of them, but first let me start with some predicting game:

## A fun problem

Take 10 coins and put them in a pile. Divide the pile into two piles of arbitrary sizes and multiply the number of coins in the two piles and record that number. Repeat this with each pile and each time record the number, until you can’t divide the piles any more. Now add up all the numbers you’ve recorded. What is the number you got?

Now my job is to guess what that number is. Well let me see, you divided the first pile and then one more time, and then another time, and… emm… I see, at the end you added 9 numbers. OK, so far, so good, then you added this one with that, and… mmm… OK… let’s see… and that adds up to… 45?

The problem is very simple if you look at it correctly. The catch is that it doesn’t matter how you divide the piles, you’ll always get the same number at the end. I think this is one of the first problems that got me interested in math. I read the problem as one of the monthly problems in a math magazine aimed at middle school students. The magazine was called Borhan, and was being published by Madreseh Publications in Iran. The editor of the Magazine, Parviz Shahriari, proposed a bunch of problems in each issue, and he would give hints or solutions in the next issue. If you want to find the actual problem take a look at the first issue. There it is 25 coins and the solution uses some basic algebra. Years later I remembered the problem again but I couldn’t remember the solutions. So, I ought to solve it myself. Here is what I came up with:

WARNING: SPOILER!

Draw 10 dots on a paper and divide it into two sets of points, say 7 and 3. Now draw all the lines that connects the dots in the first set to the second set. How many lines you’ve got? Yes, $7 \times 3$. Now keep doing this with each of the two sets of 3 and 7 dots. At the time that you can’t divide any more, you have drawn all possible lines between every pair of points. Adding up all the results of multiplications is just counting the number of lines. As you can see, it didn’t matter how you divided them. The question now is how many lines are there. Let’s number the dots: 1,2,3, …, 10. There are 9 lines that pass through 1. Put them away. Now there are 8 lines that pass through 2, put them away. Repeat this until the dot number 9, there is exactly one line. So, the total number of lines is 1+2+3+ … + 9 = 45. OK. It’s easy to do it for 10, but what if I started with 100 points? What is 1+2+3+…+99? Well, it is a famous problem due to Gauss. People say when he was in the school and his lazy teacher wanted to keep them busy after teaching them how to add, so that he could take a nap, he gives them the problem of adding all the numbers from 1 to 100. The Gauss comes up with this idea:

1 + 2 + … + 99 + 100 = S

100 + 99 + … + 2 + 1 = S

———————

101 + 101 + … + 101 + 101 = 2S

S0, there are 100 numbers that are added together and all of them are 101. The sum is 10100. But that’s twice the sum we are looking for. So let’s divide it by 2. The sum is 5050. The fact that Gauss knew how to multiply and divide right after learning how to add is not verified by anyone yet! So, now you tell me, what do you if you did it with 25 coins?

## Project-based calculus

One of the talks was about teaching calculus with projects, by Peter Zizler from Mount Royal University. There are a lot of people talking about this and there are good resources available for launching a calculus course that is entirely/mostly based on projects. To list a few

## Linear Markov Chain Markov Fields

Another talk was on Markov Chain Markov Fields by University of Calgary’s own Deniz Sezer. The talk was a review on the paper “Markov Chain Markov Field dynamics: Models and statistics” by Xavier Guyon and Cécile Hardouin, 2002 (DOI: 10.1080/02331880213192). The idea is that when dealing with Markov processes people tend to ignore the structure, because a lot of times it’s just so complicated that you really don’t know much about them. On the other hand, if you know something about the structure of what interacts with what, then you can tell a lot more about the behaviour of the system. In this talk, the focus was on the cases that only pairs interact with each other, and hence you can look at a graph. A neighbourhood of a set is defined and some properties such as ergodicity and strong law of large numbers are explored. A couple of interesting things mentioned were the notion of conditional independence of two sets $A$ and $B$ given a set $C$. This notion is related to edges of a graph to be present or not, and disconnectedness of the graph, though it wasn’t really mentioned in the talk. Another interesting concept is a set of variables $Y_A$ is called Markovian, if $A$ and $A^c$ are conditionally independent given $\partial A$ (neighbourhood of $A$). Then a theorem says

Theorem. A set of variables is Markovian if and only if all active sets in it are simplicies.

I wondered about the cases when more than just pairs can interact with each other, and hence we need something more to talk about the structure of the system, like a simplex. Apparently some few studies are done in that direction and Deniz thinks some might be found in the book “Guyon, X. (1995). Randon Fields on a Network: Modelling, Statistics and Applications (Springer)”.

## Low-rank matrix completions in geophysics

Mauricio D. Sacchi, Professor of geophysics at University of Alberta gave this interesting talk on New and not-so-new applications of low-rank matrix and tensor completions to seismic data processing. He started by mentioning the Julia Seismic software. He explained thath the Imaging is about “where” and inversion is about “what”. Then he suggested to look at the 5-dimensional data regularization. The 5 dimensions come up as 2 dimensions from the receiver, 2 dimensions from the source, and 1 is natural to be time, but it actually is the frequency. In the Parsimonious model you want to have some sparsity/simplicity, predictability, and  rank to be small (where tensor completion methods and Hankel methods play a role). The problem is that the algorithms around require regular and dense data. The solutions is to use linear and multilinear algebraic inverse theory. But there are infinitely many solutions, hence we need to put constraints (sparsity, stability, rank etc.) to get unique solutions. For example if you have

$\begin{bmatrix} d_1 \\ 0 \\ d_3 \\ d_4 \end{bmatrix} = \begin{bmatrix} 1 \\ & 0 \\ && 1 \\ &&& 1 \end{bmatrix} \begin{bmatrix} x_1 \\ x_2 \\ x_3 \\ x_4 \end{bmatrix}$,

where $d_i$s are observed data and $x_i$s are the true data, you can assume that $x = Fc$ for some Fourier transform $F$ and a sparse $c$. Then $d = SFc$. This tells us that signals that admit a sparse representation are important.

• Levy, Walker Ulrych, and Fullagar in 1982 used linear programming to solve the problem but that method is old now.
• Richard Baraniuk in 2007 used compressed sensing, but he had to assume random sampling.
• Lustig Donoho, and Pauly in 2007 also used compressed sensing for rapid MRI imaging (sparse MRI).

The problem with all of them is that they cannot guarantee e.g. 99% recovery of data.

On the other hand there is always noise. Actually the noise increases the rank, so people use rank-reduction for de-noising the data. One of the most common ways for rank reduction is by using SVD (singular value decomposition) and keeping only the largest singular values. If you find a low-rank approximation then you will have

$M^{obs} = SM$,

$M^k = M^{obs} + (I-S) R [M^{k-1}]$.

# Sagemath for calculus

I have used several geogebra applets in my calculus class, but there’s always things that you wanna do that there’s no applet out there for, or you can’t modify them etc. So, I decided to write my own code for computer algebra applications in sage, and that I can share it with students so that they can play around with it, and modify it. It makes learning more fun. I have to look into interactive option for sage, then I can actually get rid of the animations. For now, I’m just going to post a few simple codes that get the work done:
Use Sage Cloud (requires login) or Sage Cell server (does not require login) to run the following codes

# Use Newton's method to approximate a root of a function

def Newtons_method(f,a,n):
var('x')
f(x) = f
guess = a
fprime = f.derivative(x)
g(x) = x-f(x)/fprime(x)
N = n
Sguess = [guess]
print "Initial Guess: %20.16f"%(guess)
for i in range(0,N):
Nguess = g(guess)
print "Next Guess:    %20.16f"%(Nguess)
Sguess += [Nguess]
guess = Nguess.n(digits=15)

print "An estimate of the root is %20.16f"%(Sguess[N])

# Usage:
Newtons_method(x^2-57,9,5)
########################################################################

# Use Newton's method to approximate a root of a function with animation

var('x')

f(x) = x^2 - 1

xmin = 0
xmax = 3.1
ymin = -3
ymax = 9
guess = 3
fprime = f.derivative(x)
g(x) = x-f(x)/fprime(x)
N = 5
Sguess = [guess]
print "Initial Guess: %20.16f"%(guess)
P = plot(f(x), (xmin,xmax), ymin= ymin, ymax= ymax,color = "green",thickness=3)
Frames = [P]

for i in range(0,N):
P += plot(point((guess,0),size = 20,rgbcolor=(0, 0, 1)))
Frames += [P]
P += plot(point((guess,f(guess)),size = 20,rgbcolor=(1, 0, 0)))
Frames += [P]
h = f(guess) + fprime(guess)*(x-guess)
P += plot(h, (xmin,xmax),ymin= ymin, ymax= ymax, color = "blue")
Frames += [P]
Nguess = g(guess)
print "Next Guess:    %20.16f"%(Nguess)
Sguess += [Nguess]
guess = Nguess.n(digits=15)

print "An estimate of the root is %20.16f"%(Sguess[N])

animate(Frames)
########################################################################

# Finding Taylor polynomial of a function with animation

var('x')
from sage.plot.colors import rainbow
n = 20
xmin = -10
xmax = 10

g = sin(x)

P = plot(g, (xmin,xmax), color = rainbow(n+1)[0],thickness=3)
Frames = []
for i in range(n):
h = g.taylor(x, 0, i)
P += plot(h, (xmin,xmax), color = rainbow(n+1)[i+1])
Frames += [P]

m = 0
M = 12
a = animate(Frames,xmin=m,xmax=M,ymin=-1.5,ymax=1.5)
a
########################################################################
# Plot 3D with cross sections

var('x,y,z')
f = x^2 + y^2
R = 4
sum([plot3d(f,(x,-R,R),(y,level,level+0.1),frame=False) for level in srange(-R,R,1)]+[implicit_plot3d(f==z,(x,-R,R),(y,-R,R),(z,-1,32),color='khaki',opacity=0.7,frame=False)])
########################################################################
# Plot 3D with level curves

var('x,y,z')
f = x^2 + y^2
R = 4
sum([implicit_plot3d(f==level,(x,-R,R),(y,-R,R),(z,level,level+0.1),frame=False) for level in srange(0,16,2)]+[implicit_plot3d(f==z,(x,-R,R),(y,-R,R),(z,-1,16),color='khaki',opacity=0.7,frame=False)])