# Chapter 2 - Matrices

## Chapter 2.1

How can I compute with matrices?

First, the motivation.  Several of you asked whether Smiley Guy could be transformed by one thing, and then another.  Let's see!

t = var('t') @interact(layout=[['A','B'],['auto_update']]) def _(A=matrix(RDF,[[1,0],[0,1]]),B=matrix(RDF,[[1,0],[0,1]]),auto_update=False): ID = matrix(RDF,[[1,0],[0,1]]) def makeface(M): pll=M*vector((-0.5,0.5)) plr=M*vector((-0.3,0.5)) prl=M*vector((0.3,0.5)) prr=M*vector((0.5,0.5)) left_eye=line([pll,plr])+point(pll,size=5)+point(plr,size=5) right_eye=line([prl,prr],color='green')+point(prl,size=5,color='green')+point(prr,size=5,color='green') mouth=parametric_plot(M*vector([t, -0.15*sin(2*pi*t)-0.5]), (t, -0.5, 0),color='red')+parametric_plot(M*vector([t, -0.15*sin(2*pi*t)-0.5]), (t,0,0.5),color='orange') face=parametric_plot(M*vector([cos(t),sin(t)]), (t,0,pi/2),color='black')+parametric_plot(M*vector([cos(t),sin(t)]), (t,pi/2,pi),color='lavender')+parametric_plot(M*vector([cos(t),sin(t)]), (t,pi,3*pi/2),color='cyan')+parametric_plot(M*vector([cos(t),sin(t)]),(t,3*pi/2,2*pi),color='sienna') return right_eye+left_eye+face+mouth pretty_print(html('smiley guy, then transformed by $A$, and next by $B$')) G = graphics_array([[makeface(ID),makeface(A),makeface(B*A)]]) G.show(aspect_ratio=1)

## Click to the left again to hide and once more to show the dynamic interactive window

Here is Smiley Guy with rotation built in.

t = var('t') @interact def _(theta=pi/2,phi=pi/3,auto_update=False): A=matrix(RDF,[[cos(theta),-sin(theta)],[sin(theta),cos(theta)]]) B=matrix(RDF,[[cos(phi),-sin(phi)],[sin(phi),cos(phi)]]) ID = matrix(RDF,[[1,0],[0,1]]) def makeface(M): pll=M*vector((-0.5,0.5)) plr=M*vector((-0.3,0.5)) prl=M*vector((0.3,0.5)) prr=M*vector((0.5,0.5)) left_eye=line([pll,plr])+point(pll,size=5)+point(plr,size=5) right_eye=line([prl,prr],color='green')+point(prl,size=5,color='green')+point(prr,size=5,color='green') mouth=parametric_plot(M*vector([t, -0.15*sin(2*pi*t)-0.5]), (t, -0.5, 0),color='red')+parametric_plot(M*vector([t, -0.15*sin(2*pi*t)-0.5]), (t,0,0.5),color='orange') face=parametric_plot(M*vector([cos(t),sin(t)]), (t,0,pi/2),color='black')+parametric_plot(M*vector([cos(t),sin(t)]), (t,pi/2,pi),color='lavender')+parametric_plot(M*vector([cos(t),sin(t)]), (t,pi,3*pi/2),color='cyan')+parametric_plot(M*vector([cos(t),sin(t)]),(t,3*pi/2,2*pi),color='sienna') return right_eye+left_eye+face+mouth pretty_print(html('Matrices are ${}$ and ${}$'.format(latex(A),latex(B)))) pretty_print(html('smiley guy, then transformed by $A$, and next by $B$')) G = graphics_array([[makeface(ID),makeface(A),makeface(B*A)]]) G.show(aspect_ratio=1)

## Click to the left again to hide and once more to show the dynamic interactive window

Anyhow, let's define a few matrices and a vector.

V = vector([1,2]) A = matrix([[1,0],[0,-1]]) B = matrix([[0,1],[1,0]])

How do we apply a transformation to a vector?  We multiply the matrix on the left.

A*V
 (1, -2) (1, -2)

How do we apply the next transformation?  We multiply by the next matrix, on the left again.

B * (A*V)
 (-2, 1) (-2, 1)

Well, it turns out that this is associative; we could just multiply by B times A.

B*A
 [ 0 -1] [ 1 0] [ 0 -1] [ 1 0]
(B*A)*V
 (-2, 1) (-2, 1)

We can verify that multiplication doesn't always work as with numbers.

A = matrix([[0,-1],[1,0]]) B = matrix([[1,2,3],[4,5,6]]) show([A,B])
 \newcommand{\Bold}[1]{\mathbf{#1}}\left[\left(\begin{array}{rr} 0 & -1 \\ 1 & 0 \end{array}\right), \left(\begin{array}{rrr} 1 & 2 & 3 \\ 4 & 5 & 6 \end{array}\right)\right]
show(A*B)
 \newcommand{\Bold}[1]{\mathbf{#1}}\left(\begin{array}{rrr} -4 & -5 & -6 \\ 1 & 2 & 3 \end{array}\right)
show(B*A)
 Traceback (click to the left of this block for traceback) ... TypeError: unsupported operand parent(s) for '*': 'Full MatrixSpace of 2 by 3 dense matrices over Integer Ring' and 'Full MatrixSpace of 2 by 2 dense matrices over Integer Ring' Traceback (most recent call last): File "", line 1, in File "_sage_input_14.py", line 10, in exec compile(u'open("___code___.py","w").write("# -*- coding: utf-8 -*-\\n" + _support_.preparse_worksheet_cell(base64.b64decode("c2hvdyhCKkEp"),globals())+"\\n"); execfile(os.path.abspath("___code___.py")) File "", line 1, in File "/tmp/tmpZdBPFw/___code___.py", line 2, in exec compile(u'show(B*A) File "", line 1, in File "sage/structure/element.pyx", line 2949, in sage.structure.element.Matrix.__mul__ (/usr/local/sage-6.9/src/build/cythonized/sage/structure/element.c:23250) File "sage/structure/coerce.pyx", line 1070, in sage.structure.coerce.CoercionModel_cache_maps.bin_op (/usr/local/sage-6.9/src/build/cythonized/sage/structure/coerce.c:9739) TypeError: unsupported operand parent(s) for '*': 'Full MatrixSpace of 2 by 3 dense matrices over Integer Ring' and 'Full MatrixSpace of 2 by 2 dense matrices over Integer Ring'

## Chapter 2.2

Let's see what happens with Smiley Guy and the tranpose.

t = var('t') @interact(layout=[['A','B'],['auto_update']]) def _(A=matrix(RDF,[[1,0],[0,1]]),B=matrix(RDF,[[1,0],[0,1]]),auto_update=False): ID = matrix(RDF,[[1,0],[0,1]]) def makeface(M): pll=vector((-0.5,0.5))*M plr=vector((-0.3,0.5))*M prl=vector((0.3,0.5))*M prr=vector((0.5,0.5))*M left_eye=line([pll,plr])+point(pll,size=5)+point(plr,size=5) right_eye=line([prl,prr],color='green')+point(prl,size=5,color='green')+point(prr,size=5,color='green') mouth=parametric_plot(vector([t, -0.15*sin(2*pi*t)-0.5])*M, (t, -0.5, 0),color='red')+parametric_plot(vector([t, -0.15*sin(2*pi*t)-0.5])*M, (t,0,0.5),color='orange') face=parametric_plot(vector([cos(t),sin(t)])*M, (t,0,pi/2),color='black')+parametric_plot(vector([cos(t),sin(t)])*M, (t,pi/2,pi),color='lavender')+parametric_plot(vector([cos(t),sin(t)])*M, (t,pi,3*pi/2),color='cyan')+parametric_plot(vector([cos(t),sin(t)])*M,(t,3*pi/2,2*pi),color='sienna') return right_eye+left_eye+face+mouth pretty_print(html('smiley guy, then transformed by $A$, and next by $B$')) G = graphics_array([[makeface(ID),makeface(A.transpose()),makeface(A.transpose()*B.transpose())]]) G.show(aspect_ratio=1)

## Click to the left again to hide and once more to show the dynamic interactive window

Sage interprets a vector as a row or a column, whichever is appropriate.

Also, we can undo Smiley Guy's problems.  We go back to the original setup and try to do so.

t = var('t') @interact(layout=[['A','B'],['auto_update']]) def _(A=matrix(RDF,[[1,0],[0,1]]),B=matrix(RDF,[[1,0],[0,1]]),auto_update=False): ID = matrix(RDF,[[1,0],[0,1]]) def makeface(M): pll=M*vector((-0.5,0.5)) plr=M*vector((-0.3,0.5)) prl=M*vector((0.3,0.5)) prr=M*vector((0.5,0.5)) left_eye=line([pll,plr])+point(pll,size=5)+point(plr,size=5) right_eye=line([prl,prr],color='green')+point(prl,size=5,color='green')+point(prr,size=5,color='green') mouth=parametric_plot(M*vector([t, -0.15*sin(2*pi*t)-0.5]), (t, -0.5, 0),color='red')+parametric_plot(M*vector([t, -0.15*sin(2*pi*t)-0.5]), (t,0,0.5),color='orange') face=parametric_plot(M*vector([cos(t),sin(t)]), (t,0,pi/2),color='black')+parametric_plot(M*vector([cos(t),sin(t)]), (t,pi/2,pi),color='lavender')+parametric_plot(M*vector([cos(t),sin(t)]), (t,pi,3*pi/2),color='cyan')+parametric_plot(M*vector([cos(t),sin(t)]),(t,3*pi/2,2*pi),color='sienna') return right_eye+left_eye+face+mouth pretty_print(html('smiley guy, then transformed by $A$, and next by $B$')) G = graphics_array([[makeface(ID),makeface(A),makeface(B*A)]]) G.show(aspect_ratio=1)

## Click to the left again to hide and once more to show the dynamic interactive window

Here we can see what happens in general.

t = var('t') @interact def _(A=matrix(RDF,[[1,0],[0,1]]),auto_update=False): ID = matrix(RDF,[[1,0],[0,1]]) try: B = A^-1 def makeface(M): pll=M*vector((-0.5,0.5)) plr=M*vector((-0.3,0.5)) prl=M*vector((0.3,0.5)) prr=M*vector((0.5,0.5)) left_eye=line([pll,plr])+point(pll,size=5)+point(plr,size=5) right_eye=line([prl,prr],color='green')+point(prl,size=5,color='green')+point(prr,size=5,color='green') mouth=parametric_plot(M*vector([t, -0.15*sin(2*pi*t)-0.5]), (t, -0.5, 0),color='red')+parametric_plot(M*vector([t, -0.15*sin(2*pi*t)-0.5]), (t,0,0.5),color='orange') face=parametric_plot(M*vector([cos(t),sin(t)]), (t,0,pi/2),color='black')+parametric_plot(M*vector([cos(t),sin(t)]), (t,pi/2,pi),color='lavender')+parametric_plot(M*vector([cos(t),sin(t)]), (t,pi,3*pi/2),color='cyan')+parametric_plot(M*vector([cos(t),sin(t)]),(t,3*pi/2,2*pi),color='sienna') return right_eye+left_eye+face+mouth pretty_print(html('smiley guy, then transformed by $A$, and next by $A^{-1}$, which is')) show(latex(B)) G = graphics_array([[makeface(ID),makeface(A),makeface(B*A)]]) G.show(aspect_ratio=1) except: pretty_print(html("Oops! Your matrix $%s$ has no inverse!"%latex(A)))

## Click to the left again to hide and once more to show the dynamic interactive window

If you ever want to compute an inverse, use Sage.

A = matrix([ [1,2,3] , [4,5,7] , [0,8,2] ])
A^-1
 [-23/17 10/17 -1/34] [ -4/17 1/17 5/34] [ 16/17 -4/17 -3/34] [-23/17 10/17 -1/34] [ -4/17 1/17 5/34] [ 16/17 -4/17 -3/34]

Well, make sure it has one!

B = matrix(3,[1,2,3,4,5,6,7,8,9]); B
 [1 2 3] [4 5 6] [7 8 9] [1 2 3] [4 5 6] [7 8 9]
B^-1
 Traceback (click to the left of this block for traceback) ... ZeroDivisionError: input matrix must be nonsingular Traceback (most recent call last): File "", line 1, in File "_sage_input_19.py", line 10, in exec compile(u'open("___code___.py","w").write("# -*- coding: utf-8 -*-\\n" + _support_.preparse_worksheet_cell(base64.b64decode("Ql4tMQ=="),globals())+"\\n"); execfile(os.path.abspath("___code___.py")) File "", line 1, in File "/tmp/tmpiVDEES/___code___.py", line 3, in exec compile(u'B**-_sage_const_1 File "", line 1, in File "matrix0.pyx", line 4757, in sage.matrix.matrix0.Matrix.__pow__ (sage/matrix/matrix0.c:26582) File "element.pyx", line 1776, in sage.structure.element.RingElement.__pow__ (sage/structure/element.c:14627) File "element.pyx", line 3371, in sage.structure.element.generic_power_c (sage/structure/element.c:26120) File "matrix0.pyx", line 4657, in sage.matrix.matrix0.Matrix.__invert__ (sage/matrix/matrix0.c:25938) File "matrix_rational_dense.pyx", line 650, in sage.matrix.matrix_rational_dense.Matrix_rational_dense.__invert__ (sage/matrix/matrix_rational_dense.c:9364) File "matrix_rational_dense.pyx", line 756, in sage.matrix.matrix_rational_dense.Matrix_rational_dense.__invert__main (sage/matrix/matrix_rational_dense.c:9962) File "matrix_rational_dense.pyx", line 2655, in sage.matrix.matrix_rational_dense.Matrix_rational_dense._invert_pari (sage/matrix/matrix_rational_dense.c:24339) ZeroDivisionError: input matrix must be nonsingular

Also, be careful with decimals.

A = matrix(RDF,[[1,2,3],[4,5,7],[0,8,2]])
 [1.0 2.0 3.0] [4.0 5.0 7.0] [0.0 8.0 2.0] [1.0 2.0 3.0] [4.0 5.0 7.0] [0.0 8.0 2.0]
A^-1
 [ -1.352941176470588 0.588235294117647 -0.02941176470588236] [-0.23529411764705882 0.058823529411764705 0.14705882352941177] [ 0.9411764705882353 -0.23529411764705882 -0.08823529411764705] [ -1.352941176470588 0.588235294117647 -0.02941176470588236] [-0.23529411764705882 0.058823529411764705 0.14705882352941177] [ 0.9411764705882353 -0.23529411764705882 -0.08823529411764705]
A*A^-1
 [ 1.0 0.0 5.551115123125783e-17] [ 0.0 0.9999999999999996 1.1102230246251565e-16] [ 0.0 0.0 1.0] [ 1.0 0.0 5.551115123125783e-17] [ 0.0 0.9999999999999996 1.1102230246251565e-16] [ 0.0 0.0 1.0]

What the heck????  So ... is this an inverse or not?  What significance do those tiny numbers have?

This is why you need to understand the algorithms and theorems - so that things like this don't happen to you and you don't know what's going on.

## Chapter 2.3

This chapter is primarily theoretical, but there are still some interesting numerical computation issues involved.

B = matrix(RDF,3,[1..3^2]) show(B) show(B^-1)
B = matrix(ZZ,3,[1..3^2]) B^-1
 Traceback (click to the left of this block for traceback) ... ZeroDivisionError: Matrix is singular Traceback (most recent call last): File "", line 1, in File "_sage_input_10.py", line 10, in exec compile(u'open("___code___.py","w").write("# -*- coding: utf-8 -*-\\n" + _support_.preparse_worksheet_cell(base64.b64decode("QiA9IG1hdHJpeChaWiwzLFsxLi4zXjJdKQpCXi0x"),globals())+"\\n"); execfile(os.path.abspath("___code___.py")) File "", line 1, in File "/tmp/tmp68BKcW/___code___.py", line 4, in exec compile(u'B**-_sage_const_1 File "", line 1, in File "sage/matrix/matrix_integer_dense.pyx", line 1042, in sage.matrix.matrix_integer_dense.Matrix_integer_dense.__pow__ (/usr/local/sage-6.9/src/build/cythonized/sage/matrix/matrix_integer_dense.c:11506) File "sage/matrix/matrix_integer_dense.pyx", line 3851, in sage.matrix.matrix_integer_dense.Matrix_integer_dense.__invert__ (/usr/local/sage-6.9/src/build/cythonized/sage/matrix/matrix_integer_dense.c:33250) File "sage/matrix/matrix_integer_dense.pyx", line 3812, in sage.matrix.matrix_integer_dense.Matrix_integer_dense._invert_flint (/usr/local/sage-6.9/src/build/cythonized/sage/matrix/matrix_integer_dense.c:33084) ZeroDivisionError: Matrix is singular

You can get random matrices ...

random_matrix(ZZ,2,2)
 [ 0 -3] [ 0 -2] [ 0 -3] [ 0 -2]

But note that as they get bigger, there is a reason why we stick to "computer" approximations and not exact numbers.

n = 20 D = random_matrix(ZZ,n,n) B = D.change_ring(RDF) timeit("B^-1") timeit("D^-1")
 625 loops, best of 3: 71.5 µs per loop 625 loops, best of 3: 1.2 ms per loop 625 loops, best of 3: 71.5 µs per loop 625 loops, best of 3: 1.2 ms per loop

This is for # 42-44 in Section 2.3

@interact def _(x = vector([1., 2., 3., 4., 5.])): M = matrix(RDF, 5, [5, 3, 1, 7, 9, 6, 4, 2, 8, -8, 7, 5, 3, 10, 9, 9, 6, 4, -9, -5, 8, 5, 2, 11, 4]) show(M) pretty_print(html("Condition number: "+str(M.condition()))) # x = vector([1,2,3,4]) b = M*x show(b) pretty_print(html("Possible solution:")) show(M\b) pretty_print(html("Your original vector:")) show(x)

## Click to the left again to hide and once more to show the dynamic interactive window

This is a fun pattern.

@interact def _(n=[1..10]): M = matrix(n,[1 if i >= j else 0 for i in [1..n] for j in [1..n]]) show(M) show(M^-1)

## Chapter 2.4

Block matrices are implemented in Sage.

block_matrix?
 File: /usr/local/sage-5.6-linux-64bit-ubuntu_8.04.4_lts-x86_64-Linux/local/lib/python2.7/site-packages/sage/matrix/constructor.py Type: Definition: block_matrix(*args, **kwds) Docstring: Returns a larger matrix made by concatenating submatrices (rows first, then columns). For example, the matrix [ A B ] [ C D ] is made up of submatrices A, B, C, and D. INPUT: The block_matrix command takes a list of submatrices to add as blocks, optionally preceded by a ring and the number of block rows and block columns, and returns a matrix. The submatrices can be specified as a list of matrices (using nrows and ncols to determine their layout), or a list of lists of matrices, where each list forms a row. ring - the base ring nrows - the number of block rows ncols - the number of block cols sub_matrices - matrices (see below for syntax) subdivide - boolean, whether or not to add subdivision information to the matrix sparse - boolean, whether to make the resulting matrix sparse EXAMPLES: sage: A = matrix(QQ, 2, 2, [3,9,6,10]) sage: block_matrix([ [A, -A], [~A, 100*A] ]) [ 3 9| -3 -9] [ 6 10| -6 -10] [-----------+-----------] [-5/12 3/8| 300 900] [ 1/4 -1/8| 600 1000]  If the number of submatrices in each row is the same, you can specify the submatrices as a single list too: sage: block_matrix(2, 2, [ A, A, A, A ]) [ 3 9| 3 9] [ 6 10| 6 10] [-----+-----] [ 3 9| 3 9] [ 6 10| 6 10]  One can use constant entries: sage: block_matrix([ [1, A], [0, 1] ]) [ 1 0| 3 9] [ 0 1| 6 10] [-----+-----] [ 0 0| 1 0] [ 0 0| 0 1]  A zero entry may represent any square or non-square zero matrix: sage: B = matrix(QQ, 1, 1, [ 1 ] ) sage: C = matrix(QQ, 2, 2, [ 2, 3, 4, 5 ] ) sage: block_matrix([ [B, 0], [0, C] ]) [1|0 0] [-+---] [0|2 3] [0|4 5]  One can specify the number of rows or columns as keywords too: sage: block_matrix([A, -A, ~A, 100*A], ncols=4) [ 3 9| -3 -9|-5/12 3/8| 300 900] [ 6 10| -6 -10| 1/4 -1/8| 600 1000] sage: block_matrix([A, -A, ~A, 100*A], nrows=1) [ 3 9| -3 -9|-5/12 3/8| 300 900] [ 6 10| -6 -10| 1/4 -1/8| 600 1000]  It handles base rings nicely too: sage: R. = ZZ['x'] sage: block_matrix(2, 2, [1/2, A, 0, x-1]) [ 1/2 0| 3 9] [ 0 1/2| 6 10] [-----------+-----------] [ 0 0|x - 1 0] [ 0 0| 0 x - 1] sage: block_matrix(2, 2, [1/2, A, 0, x-1]).parent() Full MatrixSpace of 4 by 4 dense matrices over Univariate Polynomial Ring in x over Rational Field  Subdivisions are optional. If they are disabled, the columns need not line up: sage: B = matrix(QQ, 2, 3, range(6)) sage: block_matrix([ [~A, B], [B, ~A] ], subdivide=False) [-5/12 3/8 0 1 2] [ 1/4 -1/8 3 4 5] [ 0 1 2 -5/12 3/8] [ 3 4 5 1/4 -1/8]  Without subdivisions it also deduces dimensions for scalars if possible: sage: C = matrix(ZZ, 1, 2, range(2)) sage: block_matrix([ [ C, 0 ], [ 3, 4 ], [ 5, 6, C ] ], subdivide=False ) [0 1 0 0] [3 0 4 0] [0 3 0 4] [5 6 0 1]  If all submatrices are sparse (unless there are none at all), the result will be a sparse matrix. Otherwise it will be dense by default. The sparse keyword can be used to override this: sage: A = Matrix(ZZ, 2, 2, [0, 1, 0, 0], sparse=True) sage: block_matrix([ [ A ], [ A ] ]).parent() Full MatrixSpace of 4 by 2 sparse matrices over Integer Ring sage: block_matrix([ [ A ], [ A ] ], sparse=False).parent() Full MatrixSpace of 4 by 2 dense matrices over Integer Ring  Consecutive zero submatrices are consolidated. sage: B = matrix(2, range(4)) sage: C = matrix(2, 8, range(16)) sage: block_matrix(2, [[B,0,0,B],[C]], subdivide=False) [ 0 1 0 0 0 0 0 1] [ 2 3 0 0 0 0 2 3] [ 0 1 2 3 4 5 6 7] [ 8 9 10 11 12 13 14 15]  Ambiguity is not tolerated. sage: B = matrix(2, range(4)) sage: C = matrix(2, 8, range(16)) sage: block_matrix(2, [[B,0,B,0],[C]], subdivide=False) Traceback (click to the left of this block for traceback) ...  File: /usr/local/sage-5.6-linux-64bit-ubuntu_8.04.4_lts-x86_64-Linux/local/lib/python2.7/site-packages/sage/matrix/constructor.py Type: Definition: block_matrix(*args, **kwds) Docstring: Returns a larger matrix made by concatenating submatrices (rows first, then columns). For example, the matrix [ A B ] [ C D ] is made up of submatrices A, B, C, and D. INPUT: The block_matrix command takes a list of submatrices to add as blocks, optionally preceded by a ring and the number of block rows and block columns, and returns a matrix. The submatrices can be specified as a list of matrices (using nrows and ncols to determine their layout), or a list of lists of matrices, where each list forms a row. ring - the base ring nrows - the number of block rows ncols - the number of block cols sub_matrices - matrices (see below for syntax) subdivide - boolean, whether or not to add subdivision information to the matrix sparse - boolean, whether to make the resulting matrix sparse EXAMPLES: sage: A = matrix(QQ, 2, 2, [3,9,6,10]) sage: block_matrix([ [A, -A], [~A, 100*A] ]) [ 3 9| -3 -9] [ 6 10| -6 -10] [-----------+-----------] [-5/12 3/8| 300 900] [ 1/4 -1/8| 600 1000]  If the number of submatrices in each row is the same, you can specify the submatrices as a single list too: sage: block_matrix(2, 2, [ A, A, A, A ]) [ 3 9| 3 9] [ 6 10| 6 10] [-----+-----] [ 3 9| 3 9] [ 6 10| 6 10]  One can use constant entries: sage: block_matrix([ [1, A], [0, 1] ]) [ 1 0| 3 9] [ 0 1| 6 10] [-----+-----] [ 0 0| 1 0] [ 0 0| 0 1]  A zero entry may represent any square or non-square zero matrix: sage: B = matrix(QQ, 1, 1, [ 1 ] ) sage: C = matrix(QQ, 2, 2, [ 2, 3, 4, 5 ] ) sage: block_matrix([ [B, 0], [0, C] ]) [1|0 0] [-+---] [0|2 3] [0|4 5]  One can specify the number of rows or columns as keywords too: sage: block_matrix([A, -A, ~A, 100*A], ncols=4) [ 3 9| -3 -9|-5/12 3/8| 300 900] [ 6 10| -6 -10| 1/4 -1/8| 600 1000] sage: block_matrix([A, -A, ~A, 100*A], nrows=1) [ 3 9| -3 -9|-5/12 3/8| 300 900] [ 6 10| -6 -10| 1/4 -1/8| 600 1000]  It handles base rings nicely too: sage: R. = ZZ['x'] sage: block_matrix(2, 2, [1/2, A, 0, x-1]) [ 1/2 0| 3 9] [ 0 1/2| 6 10] [-----------+-----------] [ 0 0|x - 1 0] [ 0 0| 0 x - 1] sage: block_matrix(2, 2, [1/2, A, 0, x-1]).parent() Full MatrixSpace of 4 by 4 dense matrices over Univariate Polynomial Ring in x over Rational Field  Subdivisions are optional. If they are disabled, the columns need not line up: sage: B = matrix(QQ, 2, 3, range(6)) sage: block_matrix([ [~A, B], [B, ~A] ], subdivide=False) [-5/12 3/8 0 1 2] [ 1/4 -1/8 3 4 5] [ 0 1 2 -5/12 3/8] [ 3 4 5 1/4 -1/8]  Without subdivisions it also deduces dimensions for scalars if possible: sage: C = matrix(ZZ, 1, 2, range(2)) sage: block_matrix([ [ C, 0 ], [ 3, 4 ], [ 5, 6, C ] ], subdivide=False ) [0 1 0 0] [3 0 4 0] [0 3 0 4] [5 6 0 1]  If all submatrices are sparse (unless there are none at all), the result will be a sparse matrix. Otherwise it will be dense by default. The sparse keyword can be used to override this: sage: A = Matrix(ZZ, 2, 2, [0, 1, 0, 0], sparse=True) sage: block_matrix([ [ A ], [ A ] ]).parent() Full MatrixSpace of 4 by 2 sparse matrices over Integer Ring sage: block_matrix([ [ A ], [ A ] ], sparse=False).parent() Full MatrixSpace of 4 by 2 dense matrices over Integer Ring  Consecutive zero submatrices are consolidated. sage: B = matrix(2, range(4)) sage: C = matrix(2, 8, range(16)) sage: block_matrix(2, [[B,0,0,B],[C]], subdivide=False) [ 0 1 0 0 0 0 0 1] [ 2 3 0 0 0 0 2 3] [ 0 1 2 3 4 5 6 7] [ 8 9 10 11 12 13 14 15]  Ambiguity is not tolerated. sage: B = matrix(2, range(4)) sage: C = matrix(2, 8, range(16)) sage: block_matrix(2, [[B,0,B,0],[C]], subdivide=False) Traceback (most recent call last): ... ValueError: insufficient information to determine submatrix widths  Historically, giving only a flat list of submatrices, whose number was a perfect square, would create a new matrix by laying out the submatrices in a square grid. This behavior is now deprecated. sage: A = matrix(2, 3, range(6)) sage: B = matrix(3, 3, range(9)) sage: block_matrix([A, A, B, B]) doctest:...: DeprecationWarning: invocation of block_matrix with just a list whose length is a perfect square is deprecated. See the documentation for details. [0 1 2|0 1 2] [3 4 5|3 4 5] [-----+-----] [0 1 2|0 1 2] [3 4 5|3 4 5] [6 7 8|6 7 8]  Historically, a flat list of matrices whose number is not a perfect square, with no specification of the number of rows or columns, would raise an error. This behavior continues, but could be removed when the deprecation above is completed. sage: A = matrix(2, 3, range(6)) sage: B = matrix(3, 3, range(9)) sage: block_matrix([A, A, A, B, B, B]) Traceback (most recent call last): ... ValueError: must specify nrows or ncols for non-square block matrix. 
A = matrix([[1,2],[4,2]]) B = block_matrix([ [1, A], [0, 1] ]) B
 [1 0|1 2] [0 1|4 2] [---+---] [0 0|1 0] [0 0|0 1] [1 0|1 2] [0 1|4 2] [---+---] [0 0|1 0] [0 0|0 1]

## Chapter 2.5

Sage has LU decomposition (and many others).

M = matrix([[1,2],[3,4]]); M
 [1 2] [3 4] [1 2] [3 4]
M.LU()
 ( [0 1] [ 1 0] [ 3 4] [1 0], [1/3 1], [ 0 2/3] ) ( [0 1] [ 1 0] [ 3 4] [1 0], [1/3 1], [ 0 2/3] )

Notice that we get three matrices.  The second and third ones are the L and U matrices; the first one is one that the book refers to about possibly having to permute the rows.

show(M.LU())
 \newcommand{\Bold}[1]{\mathbf{#1}}\left(\left(\begin{array}{rr} 0 & 1 \\ 1 & 0 \end{array}\right), \left(\begin{array}{rr} 1 & 0 \\ \frac{1}{3} & 1 \end{array}\right), \left(\begin{array}{rr} 3 & 4 \\ 0 & \frac{2}{3} \end{array}\right)\right)
a,b,c = M.LU()

Here's why we need the permutation.  I hope this is clearer than the book's murky allusions.

b*c
 [3 4] [1 2] [3 4] [1 2]
a*b*c
 [1 2] [3 4] [1 2] [3 4]

Why didn't Sage just do what we did?

B = matrix([[1,0],[3,1]]) C = matrix([[1,2],[0,-2]])
show(B);show(C); show(B*C)

The reason is because in numerical routines, we would rather first switch the rows, so that the pivot is 3, not 1 (the biggest element in the column - partial pivoting).  But for our case, we can do it by hand, so we didn't rely on this technique.

LU isn't implemented for all matrices, perhaps since the underlying machinery is not available in arbitrary precision.

N = matrix([[1.,2.],[3.,4.]]); N.LU()
 Traceback (click to the left of this block for traceback) ... TypeError: base ring of the matrix must be exact, not Real Field with 53 bits of precision Traceback (most recent call last): File "", line 1, in File "_sage_input_19.py", line 10, in exec compile(u'open("___code___.py","w").write("# -*- coding: utf-8 -*-\\n" + _support_.preparse_worksheet_cell(base64.b64decode("TiA9IG1hdHJpeChbWzEuLDIuXSxbMy4sNC5dXSk7IE4uTFUoKQ=="),globals())+"\\n"); execfile(os.path.abspath("___code___.py")) File "", line 1, in File "/tmp/tmp7SemNq/___code___.py", line 3, in exec compile(u'N = matrix([[_sage_const_1p ,_sage_const_2p ],[_sage_const_3p ,_sage_const_4p ]]); N.LU() File "", line 1, in File "matrix2.pyx", line 10302, in sage.matrix.matrix2.Matrix.LU (sage/matrix/matrix2.c:49573) TypeError: base ring of the matrix must be exact, not Real Field with 53 bits of precision

But with the standard decimals in most computers, we do have it.  This uses the core LAPACK library (via Scipy) that we have mentioned before.

N = matrix(RDF,[[1.,2.],[3.,4.]]); show(N.LU())
 \newcommand{\Bold}[1]{\mathbf{#1}}\left(\left(\begin{array}{rr} 0.0 & 1.0 \\ 1.0 & 0.0 \end{array}\right), \left(\begin{array}{rr} 1.0 & 0.0 \\ 0.333333333333 & 1.0 \end{array}\right), \left(\begin{array}{rr} 3.0 & 4.0 \\ 0.0 & 0.666666666667 \end{array}\right)\right)

Let's look at Exercise 31.  I'll use exact matrices, not approximations, though in practice this would be done with decimals.

A = matrix([[4, -1, -1, 0,0,0,0,0], [-1,4,0,-1,0,0,0,0], [-1,0,4,-1,-1,0,0,0], [0,-1,-1,4,0,-1,0,0], [0,0,-1,0,4,-1,-1,0], [0,0,0,-1,-1,4,0,-1], [0,0,0,0,-1,0,4,-1], [0,0,0,0,0,-1,-1,4]]); show(A)
 \newcommand{\Bold}[1]{\mathbf{#1}}\left(\begin{array}{rrrrrrrr} 4 & -1 & -1 & 0 & 0 & 0 & 0 & 0 \\ -1 & 4 & 0 & -1 & 0 & 0 & 0 & 0 \\ -1 & 0 & 4 & -1 & -1 & 0 & 0 & 0 \\ 0 & -1 & -1 & 4 & 0 & -1 & 0 & 0 \\ 0 & 0 & -1 & 0 & 4 & -1 & -1 & 0 \\ 0 & 0 & 0 & -1 & -1 & 4 & 0 & -1 \\ 0 & 0 & 0 & 0 & -1 & 0 & 4 & -1 \\ 0 & 0 & 0 & 0 & 0 & -1 & -1 & 4 \end{array}\right)

In this case, there is no permutation needed, so we get the LU decomposition without further ado.

a,L,U = A.LU(); show(L); show(U)
b = vector([5,15,0,10,0,10,20,30])
x = U\(L\b); x
 (827/209, 1377/209, 886/209, 1546/209, 1171/209, 1831/209, 1967/209, 2517/209) (827/209, 1377/209, 886/209, 1546/209, 1171/209, 1831/209, 1967/209, 2517/209)
A*x
 (5, 15, 0, 10, 0, 10, 20, 30) (5, 15, 0, 10, 0, 10, 20, 30)

Yippee, it worked!

Notice how ugly the inverse is.  And imagine this for a much more realistic situation - not eight nodes, but tens of thousands, perhaps while you are trying to model heat changes in an airplane battery to make sure it doesn't catch fire...

show(A^-1)
 \newcommand{\Bold}[1]{\mathbf{#1}}\left(\begin{array}{rrrrrrrr} \frac{8948}{30305} & \frac{2623}{30305} & \frac{2864}{30305} & \frac{1544}{30305} & \frac{964}{30305} & \frac{689}{30305} & \frac{303}{30305} & \frac{248}{30305} \\ \frac{2623}{30305} & \frac{8948}{30305} & \frac{1544}{30305} & \frac{2864}{30305} & \frac{689}{30305} & \frac{964}{30305} & \frac{248}{30305} & \frac{303}{30305} \\ \frac{2864}{30305} & \frac{1544}{30305} & \frac{9912}{30305} & \frac{3312}{30305} & \frac{3167}{30305} & \frac{1792}{30305} & \frac{964}{30305} & \frac{689}{30305} \\ \frac{1544}{30305} & \frac{2864}{30305} & \frac{3312}{30305} & \frac{9912}{30305} & \frac{1792}{30305} & \frac{3167}{30305} & \frac{689}{30305} & \frac{964}{30305} \\ \frac{964}{30305} & \frac{689}{30305} & \frac{3167}{30305} & \frac{1792}{30305} & \frac{9912}{30305} & \frac{3312}{30305} & \frac{2864}{30305} & \frac{1544}{30305} \\ \frac{689}{30305} & \frac{964}{30305} & \frac{1792}{30305} & \frac{3167}{30305} & \frac{3312}{30305} & \frac{9912}{30305} & \frac{1544}{30305} & \frac{2864}{30305} \\ \frac{303}{30305} & \frac{248}{30305} & \frac{964}{30305} & \frac{689}{30305} & \frac{2864}{30305} & \frac{1544}{30305} & \frac{8948}{30305} & \frac{2623}{30305} \\ \frac{248}{30305} & \frac{303}{30305} & \frac{689}{30305} & \frac{964}{30305} & \frac{1544}{30305} & \frac{2864}{30305} & \frac{2623}{30305} & \frac{8948}{30305} \end{array}\right)

## Chapter 2.6

Just in case you don't believe that the series works... this is the example from #5, which we did in class.

I = identity_matrix(2) C = matrix([[0.,.5],[.6,.2]]) show((I-C)^-1)
 \newcommand{\Bold}[1]{\mathbf{#1}}\left(\begin{array}{rr} 1.60000000000000 & 1.00000000000000 \\ 1.20000000000000 & 2.00000000000000 \end{array}\right)
show(sum(C^n for n in [0..10]))
 \newcommand{\Bold}[1]{\mathbf{#1}}\left(\begin{array}{rr} 1.58835092800000 & 0.987112896000000 \\ 1.18453547520000 & 1.98319608640000 \end{array}\right)

## Chapter 2.7

Computer graphics are everywhere.  How did we make Smiley Guy?

t = var('t') @interact() def _(A=matrix(RDF,[[1,0],[0,1]])): def makeface(M): pll=M*vector((-0.5,0.5)) plr=M*vector((-0.3,0.5)) prl=M*vector((0.3,0.5)) prr=M*vector((0.5,0.5)) mouth=parametric_plot(M*vector([t, -0.15*sin(2*pi*t)-0.5]), (t, -0.5, 0),color='red')+parametric_plot(M*vector([t, -0.15*sin(2*pi*t)-0.5]), (t,0,0.5),color='orange') face=parametric_plot(M*vector([cos(t),sin(t)]), (t,0,pi/2),color='black')+parametric_plot(M*vector([cos(t),sin(t)]), (t,pi/2,pi),color='lavender')+parametric_plot(M*vector([cos(t),sin(t)]), (t,pi,3*pi/2),color='cyan')+parametric_plot(M*vector([cos(t),sin(t)]),(t,3*pi/2,2*pi),color='sienna') left_eye=line([pll,plr])+point(pll,size=5)+point(plr,size=5) right_eye=line([prl,prr],color='green')+point(prl,size=5,color='green')+point(prr,size=5,color='green') return right_eye+left_eye+face+mouth pretty_print(html('smiley guy, transformed by $A$')) G = makeface(A) G.show(aspect_ratio=1,figsize=4)

## Click to the left again to hide and once more to show the dynamic interactive window

Note how nearly everything involves a vector times a matrix.  We are moving points around with a vector.

Let's try an easier example.

@interact() def _(A=matrix(RDF,[[1,0],[0,1]])): def maketriangle(M): vertex1 = M*vector((0,0)) vertex2 = M*vector((1,3)) vertex3 = M*vector((-1,2)) edges = line([vertex1,vertex2,vertex3,vertex1])+point(vertex1,size=40,color='black')+point(vertex2,size=40,color='red')+point(vertex3,size=40,color='green')+point((0,0),size=0) return edges pretty_print(html('triangle, transformed by $%s$'%latex(A))) G = maketriangle(A) + maketriangle(identity_matrix(2)) G.show(aspect_ratio=1,figsize=4)

## Click to the left again to hide and once more to show the dynamic interactive window

We specify the vertices, and then we move them with a transformation.  THEN, and only then, do we connect the lines; this way we don't have to think about how to move the lines, just the vertices.

Problem: why can't we move it away from the origin?

Answer: we can - using block matrices and homogeneous coordinates!

@interact() def _(v = matrix([[0],[0]])): A=matrix(RDF,[[1,0],[0,1]]) def maketriangle(M,tr): M = block_matrix([[M,tr],[0,identity_matrix(1)]]) vertex1 = vector((M*vector((0,0,1)) )[:2] ) vertex2 = vector((M*vector((1,3,1)) )[:2] ) vertex3 = vector((M*vector((-1,2,1)) )[:2] ) edges = line([vertex1,vertex2,vertex3,vertex1])+point(vertex1,size=40,color='black')+point(vertex2,size=40,color='red')+point(vertex3,size=40,color='green')+point((0,0),size=0) return edges pretty_print(html('triangle, translated by $%s$'%latex(v))) G = maketriangle(A,v) + maketriangle(identity_matrix(2),matrix([[0],[0]])) G.show(aspect_ratio=1,figsize=4)

## Click to the left again to hide and once more to show the dynamic interactive window

We can combine them.  Recall that the order of translation and rotation matters!  I have created a block matrix in the code here so that we can do it just like in class.

@interact def _(v = matrix([[0],[0]]) , A = matrix(RDF,[[1,0],[0,1]]) ): def maketriangle(M,tr): M = block_matrix([[M,tr],[0,identity_matrix(1)]]) vertex1 = vector((M*vector((0,0,1)) )[:2] ) vertex2 = vector((M*vector((1,3,1)) )[:2] ) vertex3 = vector((M*vector((-1,2,1)) )[:2] ) edges = line([vertex1,vertex2,vertex3,vertex1])+point(vertex1,size=40,color='black')+point(vertex2,size=40,color='red')+point(vertex3,size=40,color='green')+point((0,0),size=0) return edges pretty_print(html('triangle, translated by $%s$ then transformed by $%s$'%(latex(v),latex(A)))) pretty_print(html('then vice versa')) G = maketriangle(A,v) + maketriangle(identity_matrix(2),matrix([[0],[0]])) H = maketriangle(A,A*v) + maketriangle(identity_matrix(2),matrix([[0],[0]])) show(graphics_array([[H,G]]),aspect_ratio=1,figsize=4)