A Programmer’s Intuition for Matrix Multiplication
betterexplained.com
betterexplained.com
Recognizing though that method of understanding something is pretty personal, I will say that what really helped me refresh was the 3Blue1Brown Essence of Linear Algebra series on Youtube. If you're trying to better grasp the subject, do yourself a favor and check out those videos.
I especially enjoyed this video (not exactly a math tutorial more of a cool explanation): https://www.youtube.com/watch?v=OkmNXy7er84
In particular his recent "Lockdown Math" series taught me a lot.
I knew and used trigonometric functions and exponentials before but they never really "clicked".
Grant is a tremendous educator.
https://youtu.be/ZxYOEwM6Wbk?t=481
Is this really the case? I checked Wikipedia and MathWorld and nobody makes a distinction between e^x and exp(x).
Even if it's a "white lie" for didactic reasons, I don't buy it, it will be much more confusing down the line for students.
Math is about finding structures governed by some rules and then generalizing them. To me the e^x defined as repeated multiplication conceptually is the same thing as the exponential function, just over broader domain. Why? You can interpolate between integer-valued X's using geometric mean. e^3 = sqrt(e^2 * e^4). What stops you from interpolating this recursively to achieve in the limit the exponential function over real numbers?
Of course, you can also go another route: you can define the number e (e.g. as the limit of (1+1/n)^n), then you can straightforwardly define natural number powers of e by repeated multiplication, then integer and rational powers through reciprocals and roots (you have to prove n-th roots exist, but that's doable), and then you can define real number exponents via limits (maybe this is what you mean by "interpolating recursively"). Now, when we come to complex numbers, you can use Euler's formula as a definition instead of a theorem, and define exp(a+bi):=e^a * (cos b + i sin b). Of course, at the end of this whole exercise you can prove that this definition is exactly equal to the power series definition.
Another benefit of the power series definition is that it also generalises e.g. to the matrix exponential (exp(A), when A is a matrix).
What I might find a little more difficult to teach to people through this lens is the phenomenon of eigenvectors, but I suppose that's to be expected—nothing will work well for all purposes.
Why? Eigenvectors are simply inputs to the network where the output keeps its shape, that is, at most it gets rescaled, as if you had applied a uniform gain to the components, but otherwise it will be the same as the input.
- does a matrix have the same left-eigenvectors as its right-eigenvectors?
- what is the relationship between the left-eigenvalues and right-eigenvalues?
- is there always an eigenvector? when is there a complete set? how do you generalize your notion of eigenvectors so that matrices always have a complete set of them?
I am not sure any of these are intuitive here.
Notice that u and v cannot be equal, because they are not the same shape. However, if v is a right eigenvector with eigenvalue \lambda, then v^T is a left eigenvector with eigenvalue \lambda, as well. More or less what this means is that we tend to just ignore left eigenvalues and only use the right eigenvalues, because the math is exactly the same up to a transpose operation.
A matrix does not necessarily have nontrivial eigenvectors. Think about the 0 matrix here. But, if A is nonzero, then it must have at least one nonzero eigenvalue, hence one nontrivial eigenvector. This is because a nonzero matrix must have at least one nonzero row and column; however, if you construct the other rows (columns) to be multiples of the first column, you end up with a matrix with only one nontrivial eigenvalue. This also illustrates the conditions necessary for an nxn matrix A to have n linearly independent eigenvectors: the rows of A must be linearly independent.
All of this is covered in a decent undergrad linear algebra course. I would suggest either finding a video course, or getting a good book and working through it, if you want to understand these things better.
Consider the matrix
M = [ 0.50 0.50 ]
[ 0.25 0.75 ]
> A matrix does not necessarily have nontrivial eigenvectors. Think about the 0 matrix here.
The zero matrix has _all_ the nontrivial vectors as eigenvectors.
> But, if A is nonzero, then it must have at least one nonzero eigenvalue, hence one nontrivial eigenvector.
Consider
N = [ 0 1 ]
[ 0 0 ]
I don't know why you consider the eigenvalue zero to somehow not count as an eigenvalue. Very suspicious.
> This also illustrates the conditions necessary for an nxn matrix A to have n linearly independent eigenvectors: the rows of A must be linearly independent.
Consider
P = [ 1 1 ]
[ 0 1 ]
> All of this is covered in a decent undergrad linear algebra course. I would suggest either finding a video course, or getting a good book and working through it, if you want to understand these things better.
sigh ... See this is why I sometimes feel bad about commenting here. Like, this comment thread was supposed to be a celebration of this different way of thinking about linear algebra, and then I have to deal with this stuff. Like I totally don’t think you meant to come across as condescending, but given that I have taught crash courses for struggling friends on linear algebra concepts they missed in their undergrad to get through our graduate work in physics, you know, it kind of does come across that way.
For the first matrix, maybe a simpler example is M = [1 1; 0 0]. Then [1; 0] is a right eigenvector with eigenvalue 1, but [1 0] M = [1 1].
Zero can of course be a proper eigenvalue. Probably they confused it with the fact that the zero vector typically doesn't count as an eigenvector.
How about this matrix?
[[0, -1],
[[1, 0]]
A snippet of code producing a counterexample:
import numpy as np
import scipy.linalg as spla
A = np.random.randn(3, 3)
right_eigenvector = spla.eig(A)[1][:,0]
right_eigenvalue = ((A @ right_eigenvector) / right_eigenvector)[0]
potential_left_eigenvector = right_eigenvector[np.newaxis, :]
# if all components of the following are not the same, then it's not
# a left eigenvector
print((potential_left_eigenvector @ A) / potential_left_eigenvector)
# prints [[-1.1327836 -0.j -0.14850693-0.j -1.84397691+0.j]]The correct theorem is that every complex square matrix has an eigenvector. If we interpret your matrix as a complex matrix, then it does have two eigenvectors, namely [1; ±i]. Hence why I’d say it kind of “doesn’t count.”
(I also think that the counterexample exhibited "counts" in the way that it's intuitively very clear why it can't have a real eigenvector. The R^2 plane can be easily visualised, as opposed to C^2, and the matrix exhibited acts as a rotation of R^2, so of course there can't be an eigenvector. As you said, once you allow complex numbers, that's not true anymore, but at that point the visual intuition breaks down somewhat.)
This description changed that. Thank you!
OTOH, it is just vectors in a many-dimensional space. We're neigh-supernatural machines for describing vectors, that's how we throw things better than any other animal, I don't really know why we need analogies here. Why do we want x'x? It maps to a distance...
..And to be clear, I am not a purist.
Joking aside, matrices are linear maps in the context of multiplication. They can absolutely have meanings assigned outside the usual "when multiplied by a vector" type of calculations. For example the "adjacency matrix" representation of a graph is very often never multiplied by a vector even in academic algorithm descriptions. YMMV if you call it a (lookup) table in that context but the term "matrix" seems fairly well established there.
Not disputing your basic point, but I'm not sure "very often never" is entirely accurate. In the GraphBLAS API vector/matrix and matrix/vector multiplication is key for almost all algorithms that depend on breadth first search, this is because matrix/vector multiplication is a single step in a BFS which is the core operation of the GraphBLAS. The vector is typically used to hold the current frontier (and mask). See the "Shortest Path Lengths" example here:
https://github.com/michelp/pygraphblas/blob/master/pygraphbl...
(disclosure, I am the pygraphblas author).
> Purist test: aren't vectors just Nx1 matrices?
I like it. :) As an interesting point, the GraphBLAS defines both Matrix and Vector types, but internally in the SuiteSparse implementation, they are all just unified into one type and yes, vectors are Nx1 or 1xN.
Not true. There is a whole field (graph signal processing) that studies graphs in terms of their description in matrix form (adjacency, connectivity and more). For example, you can find connected components by looking at the eigenvectors of a certain matrix.
Unless you're doing spectral graph theory, which looks at properties of the eigenvalues and eigenvectors of the graph's adjacency matrix. Even then, you don't usually multiply the adjacency matrix by a vector explicitly.
For instance, if A is the adjacency matrix of a graph G, with n vertices (we say 'vertices' in the math department almost exclusively -- rarely do we say 'nodes'), and k edges, denoting the trace of A by Tr(A), we can find the number of triangles in A by computing Tr(A^3)/6. We can also tell if G is connected by computing A^n. [0]
OTOH, the spectrum doesn't tell you nearly everything you'd want to know about a graph. There are many examples of non-isomorphic graphs that are nonetheless cospectral. For instance, one of my professors from grad school wrote a landmark paper in spectral graph theory in 1973, showing that almost all trees are cospectral. [1]
---
[0]: Of course, you'd probably never actually do this, because it's horribly inefficient, but here's how it goes.
First, you prove easily by induction that A^c_{i,j} is the number of paths of length c starting at vertex i and ending at vertex j. With that in mind, we see that if there is a row of all nonzero numbers in any matrix A^c for c <= n, then G is connected. With a little more work, you can even extract all the connected components of G, by looking at submatrices of A^c.
[1]: Unsurprisingly, the paper is titled "Almost all trees are cospectral". See: https://www.semanticscholar.org/paper/Almost-all-trees-are-c...
The short answer is. Yes, but if your PL treats it that way you will wind up in trouble [0].
The long answer is: https://www.youtube.com/watch?v=C2RO34b_oPM (but seriously take the time to watch this, it's fantastic).
[0] Separately C- and often C-derived languages/frameworks like numpy also treat them as 1xN matrices, but that's a nearly tabs-vs-spaces level religious war, with actual consequences. I think a lot of people got confused when translating andrew ng's early octave-based ML coursera course (Nx1) to python (1xN), and also in the back and forth between math ML papers which use Nx1 and python.
I may be misunderstanding something, but in numpy at least a vector, a 1×N array ("row vector"), and a N×1 array ("column vector") behave differently.
vec = numpy.array([1,2,3])
assert vec.shape == (3,)
assert vec.ndim == 1
row = np.array([[1,2,3]])
assert row.shape == (1,3)
assert row.ndim == 2
col = np.array([[1], [2], [3]])
assert col.shape == (3,1)
assert col.ndim == 2
assert all(vec @ numpy.eye(3) == numpy.array([1, 2, 3]))
assert all(numpy.eye(3) @ vec == numpy.array([1, 2, 3]))
The JuliaCon talk you linked does a great job of explaining what different fields mean by "vector"; thanks for linking it.
>>> import numpy as np
>>> vec = np.array([1, 2, 3])
>>> mtx = np.eye(3)
>>> np.dot(vec, mtx)
array([1., 2., 3.])
>>> np.dot(mtx, vec)
array([1., 2., 3.])
points = np.array([[0, 0, 0], [1, 0, 0], [1, 1, 0], [0, 1, 0]])
M = np.array([[0.866, 0.5, 0], [-0.5, 0.866, 0], [0, 0, 1]]) # rotation
points @ M # works
M @ points # raises ValueError, also a different transform
However, I think you can track dimensional compatibility using xarray:
import xarray as xr
import numpy as np
points = xr.DataArray([[0, 0, 0], [1, 0, 0], [1, 1, 0], [0, 1, 0]], dims=("idx", "a"))
M = xr.DataArray([[0.154, 0.988, 0], [-0.988, 0.154, 0], [0, 0, 1]], dims=("a", "b")) # rotation
points.dot(M) # works
M.dot(points) # also works
# M's a-dim is matched to point's a-dim, so both
# results have same values regardless of
# multiplication order, just transposed.
assert np.all(points.dot(M).T == points.dot(M))
In the real world, When you start doing matrix-vector, if you aren't careful with row-major vs column-major storage you can introduce performance regressions. Let's say the architecture can do 256 bit memory move into 4xf64 mmx registers, you have to think carefully about what is going where before proceeding. Some architectures might be able to do an arbitrary stride so that those f64s can be nonadjacent in a simultaneous move. I suspect not handling this correctly is part of what holds back really good fp8, fp16 performance, and I suspect for example the TPUs struggle with this.
Of course optimizing for matrix-vector can be different for optimizing for matrix-matrix, since the column vs row memory access preference for the front matrix is exactly the opposite of the back matrix.
It's been a while since I formally studied maths, but - can you describe a situation in which this is _not_ true? (modulo transposition)
So, to take a step back, I'd say that a covector is a linear function from a vector input to a scalar output (those terms themselves being somewhat fungible, the space of "vectors" needs to probably form a module over the "scalars" which must themselves be a ring—the issue is that I would as a physicist include scalar fields but you can have two nonzero fields multiply to give a zero field, so they do not obey the field axioms).
Then I would say that an [m n]-tensor is technically a multilinear function from m covectors and n vectors to a scalar. This means that a [0 1]-tensor is really exactly a covector, but a [1 0]-tensor is not really exactly a vector. One can manifestly be constructed from any vector v (take any covector c and produce the scalar c(v)) and we generally have some criteria that I forget which ensures that this is not just an injection but also a bijection implying some inverse, so that these [1 0] "co-covectors" are indeed isomorphic to the vectors. And then of course this becomes much easier in a metric space where we identify a canonical bijection between vectors and covectors called the "metric", this gives a much more direct way by which co-covectors are isomorphic to vectors.
Hopefully that all makes sense. The situations in which this is not true are non-metric spaces where linear maps from (linear maps from vectors to scalars) to scalars are just a subset of the "apply to this vector" instances, but other linear maps are maybe also included.
Then there is a more dramatic statement which I think requires the space to be Hausdorff-paracompact or something, which is the idea that of course we have the outer product which naturally combines an [a b]-tensor and a [c d] tensor into an [a+c, b+d]-tensor ("take these vectors and covectors, give them to the first tensor, take those vectors and covectors, give them to the second tensor, then take the two scalar outputs and multiply them together), and we claim that actually an [m n]-tensor can always be realized as a finite sum of outer products of m [1 0]-tensors with n [0 1]-tensors, which allows you to do arbitrary index contractions in a geometric way.
Similarly, in a vector space there can be abstract vectors, which exist no matter of what basis you choose to write them down. (Choosing a basis = choosing a way to encode them as a 1xN or Nx1 matrix).
Often in the theory, general theorems can be proven very efficiently using these abstract vectors, but most of the time when one wants to actually do things to those vectors, picking a convenient basis and working with N coordinates is the best choice.
If you tell me storing a lookup table is matrix, I think that is sort of out of context. But then again, I could be labelled a purist here, and I'm trying really hard to not fall into the trap.
Pushing one forms to two forms via exterior differentiation, for example. The chief operation is a matrix multiplication with said adjacency matrix.
IMO, just a couple of chapters from an introductory pure mathematics book is actually the most "to the point" place to start even for people who are not huge fans of math. Just the precise definition of a vector space over a field, a linear operator, a vector basis, the derivation of matrix elements, the usual related examples and some "gymnastics" on paper can make stuff "click" and be a great start for any intended application.
In Jeremy Kun's book [1] he argues that in fact data matrices _can_ be viewed as linear transformations. See e.g. section 10.9 on the SVD
> That is, we’re saying the input is R3 and the basis vectors are people:...
> By representing the ratings this way, we’re imposing the hypothesis that the _process_ of rating movies is linear in nature. ...
To solve one of these huge linear systems Ax = b , we hardly ever want [1] to go to the expense of computing and storing an explicit matrix representation of the inverse A^{-1} . Instead, we merely need to have some function f that is able to map an input vector b to some output vector f(b)=x that happens to satisfy Ax=b .
If the matrix A has sufficiently nice properties [2] we can use an efficient iterative method such as the conjugate gradient method to map the input vector b to an output vector x by solving Ax=b. To evaluate the conjugate gradient method, we don't need to necessarily have an explicit dense or sparse representation of the matrix A, we merely need a function that encodes how A acts on a vector by matrix multiplication. E.g. some function g that maps an input vector x to the output vector g(x) = Ax.
So if we've got some appropriately defined conjugate gradient method algorithm cg(g, b), we can define the function f(b) := cg(g, b)
Calling f(b) will return x that satisfies g(x) = b, i.e. f is the inverse of g. So now we're expressing linear algebra without having any explicit matrices anywhere: but our functions g and f are of course both linear operators that encode the action of matrix multiplication by A and A^{-1} respectively.
[1] Why not? It's a lot more effort to compute A^{-1} which could be used to evaluate A^{-1}b for any given b, but often we only want to know A^{-1}b for a specific b. Also, if the matrix A encodes some local property of a physical system , then A^{-1} will encode a global property -- so it will need to be a dense matrix.
[2] Sufficiently nice properties: if it's symmetric positive definite, which arises naturally in a bunch of situations. See eg https://en.m.wikipedia.org/wiki/Conjugate_gradient_method
Here is another example - cryptography analysis. Usually password length is small, usually below 1k characters. But if you increase the password length then simple old ciphers, like Vigenere, become harder to crack. Bump the key length to over 1 million bytes and suddenly your century old methods, like index of coincidence and Kasiski examination become a problem of having huge matrices to feed to analyzer. 1 million key length suddenly generates a matrix with 256M x 256M for analyzer. And all elements in that matrix are non zero.
You can define functions on matrices just like you define them on numbers. Multiplication is a function. The determinant is a function. You can even take derivatives of these functions. The derivative of the determinant function is written in terms of the trace of the matrix. Yet another matrix function. Who would have guessed that the derivative and trace are related in this way? There's really a lot of depth here.
Pure magic.
A matrix times matrix multiplication (eg 3x3 times 3x3) is just projecting the axis of one matrix onto the other: expressing the coordinate system in terms of a different coordinate system.
The 4x3 (or 4x4) matrix is needed for translation and 3d projection for a 3d triangle onto a 2d plane (screen pixels).
For translation: it's just a hack since if you write it out it's a nice and a fast hack to simulate "movement".
Then for projection you just hack the numbers in the matrix in such a way that you shrink things further apart from the origin.
TLDR: for 3d graphics matrices are used because you can hack them to represent any kind of a transformation. The reason they are used is because graphics cards are fast at performing dot products. Plus matrix multiplication doesn't require trigonometry or division so it doesn't require the slower transcendental GPU instructions (sin, cos,...).
For other uses their interpretations is different.
I always hated the dull mathematic definition of "a matrix represents a linear map" because what the matrix actually represents is completely context dependent.
A 3x3 matrix can be viewed as something that acts on the standard unit basis vectors. The action on the standard unit basis vectors can be read out by looking at the columns. (A matrix times the ith standard basis vector is just the ith column.) Since the transform is linear, all the other points in space are carried along with the action as well. This is the "linear map" viewpoint of a rotation matrix.
I show an example here https://biro.ai/we-learned-the-wrong-way-to-matrix-multiply/ where this viewpoint lets you instantly see what a rotation matrix does by looking at the columns.
The projection and translation matrices you mentioned represent restricted sets of linear maps in some vector space, but the meaning of the restriction is also captured in other mathematical structures. See affine spaces and projective spaces.
[0] https://www.youtube.com/watch?v=fNk_zzaMoSs&list=PLZHQObOWTQ...
[1] https://ocw.mit.edu/courses/mathematics/18-06-linear-algebra...
[2] https://www.youtube.com/playlist?list=PL49CF3715CB9EF31D
There are some other interpretations, for example if you think of a linear system Ax=b as a kind of “checking machine” (does x satisfy Ax=b? Yes/no, try again), then the inverse of A can produce solutions, since x=(A inverse)b.