Clojure Linear Algebra Refresher: Eigenvalues and Eigenvectors
dragan.rocks
dragan.rocks
I have to confess that, after the initial excitement around Incanter, I haven't kept tabs on it.
Last commit to develop branch: Dec 11, 2016
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.
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...
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.
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.
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."
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.
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.)
A (very imperfect) analogy might be something like the GoF book vs Peter Norvig's essay on design patterns.
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".
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:
(require '[uncomplicate.neanderthal
[core :refer [col entry nrm2 mv scal axpy copy mm dia]]
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.
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.
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.
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.
I'm more than willing to memorize some inscrutable names in exchange for this capability.
(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.
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.
Or you might happen to hate it, but, either way, I guarantee you that you will understand it quickly.
[-4 -6]
[ 3 5]
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]
The other eigenvector comes from,
[-4 -6] + [+1 0] = [-3 -6]
[ 3 5] [ 0 +1] [ 3 6]
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)
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 orthogonal coordinates. It loses this structure in the skewed coordinates. Fortunately the physicists have a very nice notation for these cases coming from their work in relativity: it is to write vectors with both "upper indices" and "lower indices". The idea is that if we're using the basis (2, -1), (1, -1) then for each vector of the basis, we find a vector which has two properties: first, it's orthogonal to all of the other vectors of the basis; second, its dot product with the vector that it corresponds to is 1. These vectors form the "dual basis".
So the vectors perpendicular to (1, -1) have the shape (t, t) for all t, we dot (2, -1) . (t, t) = 2t - t = t, so we choose t=1 and find that the dual vector to (2, -1) is (1, 1). Similarly the vectors perpendicular to (2, -1) have shape (t, 2t), dotting this with (1, -1) gives t - 2 t = -t, setting that to 1 says that t = -1, so the dual vector to (1, -1) is (-1, -2). [These mirror the expression {a, b} = {x + y, x - 2 y} above, and for good reason.]
So the idea is that we represent every vector in two ways, with its "contravariant components" v^1, v^2 such that v^1 (2, -1) + v^2 (1, -1) is our vector, and its "covariant components" v_1, v_2 such that v_1 (1, 1) + v_2 (-1, -2) is also our vector. When we do this we discover that actually we get dot products back to a nice diagonal form even in a skewed coordinate system, that form is u_1 v^1 + u_2 v^2 + ... = u^1 v_1 + u^2 v_2 + ... . If you've got a nice crystal structure in physics, you might have, say, that the electrical conductivity really "wants" to be expressed in these skewed coordinates, where if you go down any of the natural directions of the crystal the material responds via Ohm's law. If it has slightly different resistances in slightly different crystal directions, say, then because the crystal is a skewed coordinate system, if you apply an electric field in an arbitrary physical direction you will in general create a current in a slightly different physical direction. But in the skewed coordinates it just looks like `J = {s1 E1, s2 E2, s3 E3}` corresponding to a diagonal conductivity tensor `diag(s1, s2, s3)` ... it's just that when you're not in the crystal structure you can't quite see that this material wants to flow in a couple of skewed directions preferentially.