This is nice, but I wish he showed it in mathematical syntax too. I don't understand some of the calculations because I don't know much about Clojure. Or Linear Algebra.
The idea is that you follow this with a math Linear Algebra textbook. There is a recommended textbook, but any will do. I am not a mathematician, and there is not enough space in a blog post do develop everything step by step anyway. This should be a math->programming link, not a self-contained book.
For the syntax, just take a little time to review the Clojure syntax. I guarantee you that the syntax is as simple as it gets for a programming language!! You might happen to like it.
Or you might happen to hate it, but, either way, I guarantee you that you will understand it quickly.
So the transform is what we'd write as a matrix as:
[-4 -6]
[ 3 5]
This matrix corresponds to the following linear transform: if you hand in the pair (x, y) it gives you the pair T(x, y) = (-4x - 6y, 3x + 5y). The word "linear" means that T(x1 + x2, y1 + y2) = T(x1, y1) + T(x2, y2) for all x1, x2, y1, y2 and that T(k x, k y) = k T(x, y) for all k. This has a nice interpretation in terms of the "vector space" underneath; a linear transform "distributes over vector addition" and "commutes with scalar multiplication."
In fact vector spaces are super-general. The set of infinite sequences of real numbers `[a0, a1, a2, ...]` is a very nice vector space, and the Fibonacci recurrence F[n] = F[n - 1] + F[n - 2] is linear in such infinite sequences. It can be solved by F[n] = p^n for two bases p=b1, p=b2, and then you can see that any other sequence obeying that recurrence relation must have the form F[n] = A b1^n + B b2^n for some constants A, B. It's the same principles of linear algebra which get you there. (Exercise: solve the Fibonacci recurrence for b1, b2 and find A and B such that F[0] = 0, F[1] = 1, the normal Fibonacci numbers.)
Anyway, back to the matrix. As you can see the trace of this matrix (sum of diagonal entries) is 1, and the determinant (for a 2x2 matrix, the product of the diagonal entries minus the product of the other two entries) is -20 + 18 = -2. It turns out that the trace is always the sum of the eigenvalues and the determinant is always the product of them, giving a nice way to understand this as the eigenvalues +2 and -1, the only two numbers that you can multiply together to get -2 and sum together to get +1.
If we know that these are the eigenvalues then we would calculate the eigenvectors by looking for a non-trivial nullspace, T(x) = k x implying that T2(x) = T(x) - k x maps this nontrivial vector x to 0. Since the - k x is the diagonal matrix diag(-k, -k), we are looking at for k = +2 the matrix
[-4 -6] + [-2 0] = [-6 -6]
[ 3 5] [ 0 -2] [ 3 3]
Now that looks like a very degenerate transform! T(x, y) = (-6 x - 6 y, 3 x + 3 y), what's going to map this to (0, 0)? Just x = -y will do perfectly nicely. So (1, -1) is one representative eigenvector for an entire eigendirection (t, -t) for all t. Plugging that into the original we also see T(1, -1) = (-4 + 6, 3 - 5) = (2, -2), proving that this is all 100% self-consistent.
The other eigenvector comes from,
[-4 -6] + [+1 0] = [-3 -6]
[ 3 5] [ 0 +1] [ 3 6]
This takes a little more thinking, the eigenvector is (2, -1) representing the eigendirection (2t, -t) for all t.
Now what's the basic reason that you're doing all of this? Here's what: the eigenbasis. These vectors are not orthogonal, but they do span the space with a skewed coordinate system. Any vector
(x, y) = a (2, -1) + b (1, -1)
Call these new components {a, b} with curly braces.
To solve for it explicitly, notice that {1, -1} = (1, 0) and {1, -2} = (0, 1). So x (1, 0) + y (0, 1) = x {1, -1} + y {1, -2} = {x + y, -x - 2y}, and we can find that a = x + y, b = -x -2y.
Now if you pay the pain of using these skewed coordinates to begin with, then the matrix has a very nice representation in these coordinates: T{a, b} = {-a, 2b}. It just treats both of these coordinates independently, scaling them without mixing them. If you want to repeat it n times, it will just be T^n {a, b} = {(-1)^n a, 2^n b}. (Exercise: one of those roots b1, b2 for the Fibonaccis, let's say b2, has absolute magnitude less than 1, hence its exponentiations drift towards 0. Program an algorithm to calculate the Nth Fibonacci as `round(A*pow(b1, n))` and see how it does.)
What complicates things a little is that usually your vector space is also an inner product space, and often that inner product has a nice structure (a, b) . (x, y) = a x + b y, in the orthog...
Anyone know what the status of Incanter is these days? It seems like Incanter could benefit mightily from the Neanderthal library referenced in the article.
I have to confess that, after the initial excitement around Incanter, I haven't kept tabs on it.
As far as I understand it, there's core.matrix which is intended as sort of a universal Clojure API for linear algebra. There are different implementations of the core.matrix API, like the pure JVM vectorz-clj, BLAS based clatrix and even a clojurescript one, aljabr. Incanter uses core.matrix and can therefore use these different implementations. However, the Neanderthal library does not support core.matrix, so Incanter can't use it.
As for why Neanderthal doesn't support core.matrix, I don't know. I see the author is active in this thread, so he can probably give a much better answer, but from reading the mailing lists I gather he mainly created Neanderthal to scratch his own itch, so core.matrix support wasn't a priority.
Disclaimer: my brother is mikera, maintainer of core.matrix
I asked why neanderthal wasn't part of core.matrix and Mike replied that he'd like if it were 'for the sake of the ecosystem', and he'd offered some code and tests last year to get this started but it didn't get taken in. It's still there from last October, probably a bit out of date by now...
Because when you string a couple of those please-compute-eigen-vectors-and-also-eigen-values in nested calls, you easily end up with a multi-page text in no time.
Besides, the standard way to call eigenvectors/eigenvalues function is ev in LAPACK - sgeev/dgeev.
You might have a different opinion once you have written a few tens of thousands lines of numerical code. There is a reason why mathematical notation is so terse - much terser than the “garbage” you bemoan, in fact.
Also, regarding your specific example, `eigen` would be ambiguous. Do you mean eigenvalues? Eigenvectors? Eigenfunctions? Eigenspaces?
You are walking into a classic argument about concise notation for experts versus self-explanatory names for generalist programmers. See https://news.ycombinator.com/item?id=14239228 for a previous iteration of this argument on this site.
Suffice it to say that anyone with sufficient mathematical expertise to actually make sense of this stuff, who is going to do anything interesting will want to use concise names. And yes, this will be the first thing that a programmer who doesn't know it will complain about, but it is far from the largest problem that you have to face in using this stuff. And if you work through the actually important barriers to doing useful stuff with this knowledge, you'll probably appreciate concise notation.
It's still an interesting idea to try to expose an unknown domain in terms of a well thought out API which takes care of some of those dirty details. An ORM would be an example. Should linear algebra API be the raw numerical algorithms (Fortran like) or something dressed up which throw exceptions on failed convergence?
(This is probably orthogonal to your comment btw )
Any abstraction layer built on top of linear algebra intended for non-mathematicians should be given mnemonic names that are domain specific for what that abstraction layer is aimed at.
Any abstraction layer built with the purpose of making it easy to do linear algebra in software should use names that are domain specific to linear algebra. Which will be the short concise names that make sense to those with mathematical expertise, no matter what generalist programmers think of said names.
I think most programmers probably either need control of the nitty-gritty of the numerical algorithms or would be better-served by a library explicitly modeling their domain (with library authors minding the p's and q's of the numerical algorithm).
These are not names designed for people who understand linear algebra from an undergraduate mathematics course, they’re just names taken from BLAS, which were chosen based on the conventions of Fortran code.
There’s no particularly good reason to shorten `norm` to `nrm` except to be cryptic.
Knowing that `mv` means “function that takes two arguments and multiplies the matrix in the first argument by the vector in the second argument” and that `axpb` means “takes three arguments, and multiplies a scalar by a vector and then adds it to another vector” and that `scal` means “takes two arguments and multiplies the scalar first argument by every element of the vector or matrix in the second argument” requires familiarity with BLAS. In general the names are fairly un-systematic. Once you’ve learned them reading code is reasonably straight-forward (albeit much less clear than code written using more explicit syntax like you’d find in Matlab or similar), but writing code is going to require either frequent trips to the docs or some time spent on memorization.
“Systematic” would be e.g. `ab`, `ax`, `axpy`, `kx`, instead of `mm`, `mv`, `axpy`, `scal`. (Or whatever other preferred system, hopefully something carefully thought through rather than slapped together.)
The names are obviously meaningful, they’re just fairly ad-hoc and not guessable a priori.
I would say mm and mv name the types of the arguments, like a kind of Hungarian notation, which makes them okay. I think `ab`, `ax`, `axpy` make the problem worse: in the first two a and b are matrices and in the last a is a scalar. As for `kx`, anything named `k` is usually supposed to be an integer like everywhere else in mathematics (why not `ax` then?). As for scal-scale and nrm-norm, I have no idea why they did that, it seems wrong to me too.
Furthermore, it's a specialist expert subdomain—the names are already so standard that everybody who needs them already knows them or looks them up in a second. So I don't see "not guessable" as a valid criticism.
I'm fond of "eig" which in Matlab is the function to return eigenvectors and eigenvalues. My experience is that it's easier to recognize and remember abbreviations rather than acronyms.
It would be interesting to have a language with first-class(ish) support for "tagged" aliases, where an IDE could then expand to short/long versions as configured. To some extent you could already do it in a lot of languages which support annotations/decorations, but it might be nice to have it more built into the language for "standardization" if nothing else.
What he gives is the identity matrix—multiply on the right by any 2x1 vector and you will get that vector back. So any such vector, provided it is not the zero vector, is an eigenvector with eigenvalue 1.
I obviously formulated this wrong. Will rephrase, but doesn't the intuiton generally hold that when you find an eigenvector for one eigenvalue, then you can construct infinite number of eigenvectors by scaling that one?
Yes, that is trivially true. In fact any linear combination of eigenvectors for the same eigenvalue will be an eigenvector for that eigenvalue.
Another important point to appreciate is that eigenvalues and eigenvectors frequently require complex numbers. For example the eigenvalues of [[0 1] [-1 0]] are i and -i. This geometrically correlates with some kind of "rotation".
[edit: the rest of my original comment was not correct]
If an eigenvalue is unique, the corresponding eigenvectors are on a line. If it appears twice, the corresponding eigenvectors are [edit: could be] on a plane. In general, they span [edit: could span] a subspace of dimension equal to the multiplicity.
The dimension of the space spanned by the eigenvectors associated with a degenerate eigenvalue is bounded above by the algebraic multiplicity of the eigenvalue and bounded below by one.
You seem to have reworded to "All eigenvectors that correspond to one a unique eigenvalue lie on the same line, but have different magnitudes." This is technically true if you interpret "unique" to mean "multiplicity of one", but still seems confusing, especially as you have not yet discussed the multiplicity of eigenvalues.
I would consider rephrasing the article as you do in the above comment, "when you find an eigenvector for one eigenvalue, then you can construct infinite number of eigenvectors by scaling that one".
The book based on this paper is one of the best text books I've read. I highly recommend reading it if you're looking for a solid understanding of linear algebra.
As someone who worked their way through Strang (we did a little study group and did the homework assigned to some class at Rutgers), I'd love it if someone could take a moment to explain what they like so much about Axler. I've tried a couple times to slog my way through Axler and had not much success.
I have both Axler and Strang. I learned from Axler but occasionally use Strang as a reference. I find that Axler emphasizes the algebra part of linear algebra while Strang emphasizes the linear part. One isn't better than the other but one might make more sense to different people.
Axler avoids determinants with a passion. For good reason, too, because they are not helpful in understanding how spaces interact. I forget how Strang proves the spectral theorem but compared to Lax the pedagogy could not be more different. Determinants are really nice tools for short proofs, but not helpful in understanding. But Axler is the way you want to think about linear algebra as a practitioner.
Axler is easier to digest if you have some familiarity with the tools and the jargon but hated/forgot everything in your college linear algebra course.
A (very imperfect) analogy might be something like the GoF book vs Peter Norvig's essay on design patterns.
I’m not in the target audience and have neither learned nor taught from Axler’s book, so I can’t tell you precisely what people like about it, but it is aimed at pure math undergraduates who are interested in using linear algebra in much more abstract contexts, and focuses more on proofs than many introductory linear algebra textbooks.
You don’t think this book is aimed at pure mathematics students?
It’s certainly not aimed at numerical analysis students, or engineering students, or physics students. (Which isn’t to say that those students can’t take pure math courses if they want.)
I mean in the context of Strang. I don't think Strang is aimed strictly at, say, engineering students. To me, the difference is the starting point, more than anything else.
As far as I can tell Strang’s target audience is something like: most undergraduates at MIT who didn’t already learn the subject before arrival, except the pure math students who are likely to substitute a more theoretical course (18.700 or 18.701 vs. 18.06).
Right. Each is aimed at somewhat different audiences. We're undoubtedly getting into some pretty fine hairsplitting, the thing I was whining about is 'pure maths students' and 'everyone else' is not an accurate way to describe them.
Not a contradiction. Their first exposure would have been something like multivariable Calculus. So students have seen matrices, but they were magic, unmotivated, and probably didn't make much sense except in a "here is the formula to memorize" kind of way.
Axler tries to teach students how to understand linear algebra.
If all that you want to do is use it, the prospect of that understanding may not be very motivating.
I was commenting about what Axler is, and not comparing it to anything else.
My point is that the comment that you quoted from the preface in no way changes the fact that the book's point is to convey an understanding of linear algebra that is primarily of interest to people going on in math.
Now I happen to think it is the right way to understand linear algebra and is how people in other fields should think about it. Because it is easier to figure out again if you've not done it in a while. But this point of view is primarily going to motivate would be mathematicians.
It is probably logically unnecessary but I think it's best to have experience with the "computational" aspects of linear algebra before diving into Axler's book.
I really appreciate that series, but my one criticism is that it does not connect concepts to how they are used. Maybe it's just me, but I get a much stronger intuition if I imagine what's happening when I apply the concepts, and I think it's just low hanging fruit for him to enshrine that series as possibly one of the best resources out there if he took a few minutes here and there to walk through details of application.
Yes, you're right, if we take this course as standalone.
When used alongside a more "traditional" course with applications, it's how it becomes very interesting.
You can easily feel guilty as a student when watching this, thinking you should focus on what's been specifically taught to you (during you 'official' lectures), and not on the conceptual ideas behind it. But it does make learning easier, and maybe encourages you to come back to the subject later, because you don't only see the subject as a tool anymore, and get a feeling a familiarity.
74 comments
[ 4.7 ms ] story [ 81.6 ms ] threadOr you might happen to hate it, but, either way, I guarantee you that you will understand it quickly.
In fact vector spaces are super-general. The set of infinite sequences of real numbers `[a0, a1, a2, ...]` is a very nice vector space, and the Fibonacci recurrence F[n] = F[n - 1] + F[n - 2] is linear in such infinite sequences. It can be solved by F[n] = p^n for two bases p=b1, p=b2, and then you can see that any other sequence obeying that recurrence relation must have the form F[n] = A b1^n + B b2^n for some constants A, B. It's the same principles of linear algebra which get you there. (Exercise: solve the Fibonacci recurrence for b1, b2 and find A and B such that F[0] = 0, F[1] = 1, the normal Fibonacci numbers.)
Anyway, back to the matrix. As you can see the trace of this matrix (sum of diagonal entries) is 1, and the determinant (for a 2x2 matrix, the product of the diagonal entries minus the product of the other two entries) is -20 + 18 = -2. It turns out that the trace is always the sum of the eigenvalues and the determinant is always the product of them, giving a nice way to understand this as the eigenvalues +2 and -1, the only two numbers that you can multiply together to get -2 and sum together to get +1.
If we know that these are the eigenvalues then we would calculate the eigenvectors by looking for a non-trivial nullspace, T(x) = k x implying that T2(x) = T(x) - k x maps this nontrivial vector x to 0. Since the - k x is the diagonal matrix diag(-k, -k), we are looking at for k = +2 the matrix
Now that looks like a very degenerate transform! T(x, y) = (-6 x - 6 y, 3 x + 3 y), what's going to map this to (0, 0)? Just x = -y will do perfectly nicely. So (1, -1) is one representative eigenvector for an entire eigendirection (t, -t) for all t. Plugging that into the original we also see T(1, -1) = (-4 + 6, 3 - 5) = (2, -2), proving that this is all 100% self-consistent.The other eigenvector comes from,
This takes a little more thinking, the eigenvector is (2, -1) representing the eigendirection (2t, -t) for all t.Now what's the basic reason that you're doing all of this? Here's what: the eigenbasis. These vectors are not orthogonal, but they do span the space with a skewed coordinate system. Any vector
Call these new components {a, b} with curly braces.To solve for it explicitly, notice that {1, -1} = (1, 0) and {1, -2} = (0, 1). So x (1, 0) + y (0, 1) = x {1, -1} + y {1, -2} = {x + y, -x - 2y}, and we can find that a = x + y, b = -x -2y.
Now if you pay the pain of using these skewed coordinates to begin with, then the matrix has a very nice representation in these coordinates: T{a, b} = {-a, 2b}. It just treats both of these coordinates independently, scaling them without mixing them. If you want to repeat it n times, it will just be T^n {a, b} = {(-1)^n a, 2^n b}. (Exercise: one of those roots b1, b2 for the Fibonaccis, let's say b2, has absolute magnitude less than 1, hence its exponentiations drift towards 0. Program an algorithm to calculate the Nth Fibonacci as `round(A*pow(b1, n))` and see how it does.)
What complicates things a little is that usually your vector space is also an inner product space, and often that inner product has a nice structure (a, b) . (x, y) = a x + b y, in the orthog...
I have to confess that, after the initial excitement around Incanter, I haven't kept tabs on it.
As for why Neanderthal doesn't support core.matrix, I don't know. I see the author is active in this thread, so he can probably give a much better answer, but from reading the mailing lists I gather he mainly created Neanderthal to scratch his own itch, so core.matrix support wasn't a priority.
A community member of ours made an implementation of core.matrix backed by nd4j a while back: https://github.com/ds923y/nd4clj
It needs a bit of updating though.
We support cpu and gpu cublas and arbitrary matrix operations. We have our own memory management engine among other things built in.
This lib will also support automatic differentiation and sparse in the next few months.
I asked why neanderthal wasn't part of core.matrix and Mike replied that he'd like if it were 'for the sake of the ecosystem', and he'd offered some code and tests last year to get this started but it didn't get taken in. It's still there from last October, probably a bit out of date by now...
https://github.com/mikera/neanderthal/network
Last commit to develop branch: Dec 11, 2016
https://github.com/incanter/incanter/commits/master
https://github.com/incanter/incanter/commits/develop
Why not call that function eigen-vectors instead of ev?
Some library authors need to be smacked on the head with a programming book and learn how to design sensible API's.
Seriously, Clojure is typically easy on the eyes but this is just garbage:
axpy? mm? dia? nrm2?Besides, the standard way to call eigenvectors/eigenvalues function is ev in LAPACK - sgeev/dgeev.
And even so, we know today it's preferable to have longer but easier to read code than this kind of garbage.
> sgeev/dgeev.
That's just as bad.
On the other hand, you are free to choose libraries that you use. There is surely no shortage of libraries than invent those "proper" names.
Also, regarding your specific example, `eigen` would be ambiguous. Do you mean eigenvalues? Eigenvectors? Eigenfunctions? Eigenspaces?
Those are the standard BLAS names for those things.
Suffice it to say that anyone with sufficient mathematical expertise to actually make sense of this stuff, who is going to do anything interesting will want to use concise names. And yes, this will be the first thing that a programmer who doesn't know it will complain about, but it is far from the largest problem that you have to face in using this stuff. And if you work through the actually important barriers to doing useful stuff with this knowledge, you'll probably appreciate concise notation.
(This is probably orthogonal to your comment btw )
Any abstraction layer built with the purpose of making it easy to do linear algebra in software should use names that are domain specific to linear algebra. Which will be the short concise names that make sense to those with mathematical expertise, no matter what generalist programmers think of said names.
I'm more than willing to memorize some inscrutable names in exchange for this capability.
There’s no particularly good reason to shorten `norm` to `nrm` except to be cryptic.
Knowing that `mv` means “function that takes two arguments and multiplies the matrix in the first argument by the vector in the second argument” and that `axpb` means “takes three arguments, and multiplies a scalar by a vector and then adds it to another vector” and that `scal` means “takes two arguments and multiplies the scalar first argument by every element of the vector or matrix in the second argument” requires familiarity with BLAS. In general the names are fairly un-systematic. Once you’ve learned them reading code is reasonably straight-forward (albeit much less clear than code written using more explicit syntax like you’d find in Matlab or similar), but writing code is going to require either frequent trips to the docs or some time spent on memorization.
The names are obviously meaningful, they’re just fairly ad-hoc and not guessable a priori.
xpy uses x and y because x and y could be vectors, but can also be matrices
mm can not be used for vectors, only for matrices.
Of course, the main reason is that these are names used for the last 40+ years in the standard for this type of software, BLAS.
Furthermore, it's a specialist expert subdomain—the names are already so standard that everybody who needs them already knows them or looks them up in a second. So I don't see "not guessable" as a valid criticism.
The eigenvectors for [[1 0] [0 1]] are 2 dimensional.
I highly recommend http://www.axler.net/DwD.html for developing a good intuition about what eigenvalues and eigenvectors actually are.
Another important point to appreciate is that eigenvalues and eigenvectors frequently require complex numbers. For example the eigenvalues of [[0 1] [-1 0]] are i and -i. This geometrically correlates with some kind of "rotation".
[edit: the rest of my original comment was not correct]
If an eigenvalue is unique, the corresponding eigenvectors are on a line. If it appears twice, the corresponding eigenvectors are [edit: could be] on a plane. In general, they span [edit: could span] a subspace of dimension equal to the multiplicity.
The generalized eigenvectors are on a plane.
The case where the eigenvectors do not span the space corresponds to matrices whose Jordan normal form has upper triangular bits.
I would consider rephrasing the article as you do in the above comment, "when you find an eigenvector for one eigenvalue, then you can construct infinite number of eigenvectors by scaling that one".
A (very imperfect) analogy might be something like the GoF book vs Peter Norvig's essay on design patterns.
Not sure if they’re useful (I haven’t watched them), but Axler made a series of videos about the core content of his book https://www.youtube.com/watch?v=lkx2BJcnyxk&list=PLGAnmvB9m7...
I don't think that's really true. The book itself puts it much more simply and directly, at the very beginning in 'Preface for the Instructor':
"You are about to teach a course that will probably give students their second exposure to linear algebra."
It’s certainly not aimed at numerical analysis students, or engineering students, or physics students. (Which isn’t to say that those students can’t take pure math courses if they want.)
Axler tries to teach students how to understand linear algebra.
If all that you want to do is use it, the prospect of that understanding may not be very motivating.
My point is that the comment that you quoted from the preface in no way changes the fact that the book's point is to convey an understanding of linear algebra that is primarily of interest to people going on in math.
Now I happen to think it is the right way to understand linear algebra and is how people in other fields should think about it. Because it is easier to figure out again if you've not done it in a while. But this point of view is primarily going to motivate would be mathematicians.
Sure, but the topic is 'Given Strang, what's the deal with Axler'. It's a perfectly sensible question in its own right.