the second section of course of course was on a applications in geometry in
approximations statistics areas like that. And this last section in the course is on.
Numerics and algorithms how do you actually solve these problems?
Now, you have been solving in these problems more well for the
last two parts of the course you have been solving these problems.
but presumably
at this point you don't know how it's done and the main
goal of this section, the whole section of the course is to demystify.
How it is that we solve these convex problems, right.
So that's that's the idea.
Now, most of you are not going to write your own solvers.
Right.
Actually, some of you might have to, if you solve, if
you go into areas where the problems are absolutely gigantic, you know,
something like CVX is not going to work,
you're going to have to develop your own solver.
or if you go into like real time systems where you
need to be able to solve optimization problems in micro seconds.
Again something like cvx is not going to work for you
and you're going to really have to known how its done.
I should say a couple more things about that.
I still believe strongly, its not just a question of demystifying it
but I believe its very important for everyone to know a little bit
About how these problems are solved and specifically what you want to know
is the link between problem structure and how fast we can solve it.
So that's really, if there was one take away for
the entire last part of this course, it's to understand that.
That one thing.
You should be able to look at that problem and say "Ooh!
That's got a hundred thousand variables, but I can solve that super fast, because
I recognize the structure." So that's what this whole section of the course is about.
Now we are going to go back today and start, well, at the very beginning.
We are going to talk a little about numerical algebra.
Some of you probably know this, others of you probably don't.
People who have used things like...
Matlab have certainly been using this as a user, as
a consumer but you may not know actually how it's done.
and what we'll do today is just look at the basics
and I should also say how it fits into the big picture.
So the big picture is, this lecture, we're going to
talk about How do you solve a linear equation.
So that's fine.
And there's a finite algorithm for that and we'll look at
how that's done, and we'll we'll see we'll see how it's done.
It's gaussian elimination.
You've been doing this since high school I assume.
Okay then I'll tell you how this fits into the rest of the course.
Right?
So in the rest of the course we'll see how do you
minimize, we'll start by talking how do you minimize a smooth function.
well, with start with how do you, we'll
start with how do you minimize a quadratic function.
Well, you set the gradient equal to 0.
But the grading of a quadratic fuction is a linear
function, and so, minimizing a quadratic function, solving linear equations.
So if you like, you can think of this whole stuff we're going to look at now,
how do you solve linear equations, as how do you minimize a quadratic funtion.
So, that's the optimization connection. Okay so, we will first start
by talking matrix structure and algorithm complexity.
matrix structure is, you know, somethings are obvious like
what is, you know, diagonal matrices, identity matrices sparse matrices.
We'll talk about that.
then we'll talk about the dominant method used to solve linear equations.
and that's using factored matrices.
So you'll factor a matrix and then solve, do a sequence
of solves, with the factored, with the factors.
we'll talk about some very famous factorizations, the so called
LU, that's going to be for lower triangular, upper triangular, the
Cholesky factorization and something called
the LDL transpose factorization where L
is lower triangular, D is diagonal and L is upper triangular.
Actually we will see in some cases D is block diagram.
Then, we'll look at a general construction called block elimination
in the matice inversion lemma and this is something that shows you how
it is that you can solve
sometimes absolutely enormous sets of linear equations.
very efficeintly by explointing some structure that's present.
And the last thing we'll talk about,
though very briefly is solving undetermined equations.
So first let's just start with solving Ax equals b.
where A is n by n.
Okay?
And in fact let's say on that laptop right there,
that's not mine, but I know the, I know the number.
I wa, I want to ask you all, how long do
you think it takes to solve the equation Ax equals b?
Where x is size 1,000.
Before you say anything, I want to point a few things out.
That's solving
1,000 equations, for 1,000 variables. The size of A, is a million coefficients.
Everybody got it? Okay?
And now what I want to know is how fast can I solve it?
The answer is 30 milliseconds.
Okay, and I happen to know the number on that thing is about ten.
So actually, this is an amazing thing, and it's actually
worthwhile for all of us to sit and think how unbelievably amazing that is, right.
i mean that is ridiculous.
Right?
That's something that takes like a billion, we'll
see it takes, generic methods take a billion operations.
It's ten, fifteen miliseconds on that. Okay?
It's got multiple cores, does all sort of things are happening to make that number.
But just the idea actually that it's less than a second is unbelievably impressive.
Okay? So well we'll get to that.
So it, so it basically for general methods right?
The general methods means you, you're exploiting
nothing in the matrix a, absolutely none, right?
All the entries are assumed to be non-zero there's no particular structure.
The structure would be something like a might be
A discreet foray transform or a discreet cosine transform.
In which case, of course, you can solve it faster than n cubed.
You can solve it in like n log n. Okay so that would be structure.
Other structure would be sparsity.
Right, that would be that A has a lot of zeroes.
Right, one definition of sparsity goes something like this: it's.
A matrix that has enough zeros that it's worth your while to take advantage of it.
That's 1967; I mean it's a perfectly good definition, right?
So,
okay.
So general methods, it says that these things grow, like N cubed.
oh, so for fun, how long would it take that laptop there to solve?
Say a hundred, hundred equations,ah and a hundred
variables that's 10 thousand coefficients in the matrix.
So the answer was what 30? >> Micro seconds
>> Micros seconds, thank you, okay. So I mean that's ridiculous,
by the way this also hints Why you should know about all those things, right.
These numbers are so ridiculous that it means these things can now be embedded
on all, it means the amount of computation you can do can now, is unbelievable.
you can put stuff, you can actually embed, all of these things can be embedded.
They can, and things that would normally take like seconds to run.
>> can run in milliseconds. They can be embedded in things.
I mean, crazy stuff
can happen. So, OK.
all right.
So, basically for general, these are the generic methods.
This grows as N cubed.
And I'll, I mean I'll tell you what the generic methods are, shortly.
what, what they are.
and this is, this is the part that I think.
I want everybody to know.
when A is structured, and we'll talk about all sorts of structure, it's way faster.
Or generally, substantially faster than n cubed.
Okay, and structured is stuff like banded.
A banded matrix is 1 Where, basically, a i j is
zero, if the difference between i and j is big enough.
Bigger than some thresh hold.
Right, systems like that can actually be solved in linear time, okay.
So, I mean these are just things you did, you wouldn't I mean unless you knew this,
it's not obvious right, so what that say is I can solve a million variables with
a million equations on that laptop, actually I can do that on my phone
but on that,that laptop I solve a
system of a million equations,a million variables.
If, if it's banded with a bandwidth of like five.
Right?
And not only that, it can solve, I mean, unbelievably fast.
I won't do the calculation now but it'll be shocking.
OK?
And if you didn't know that, you'd think
well that's, that's a big hard problem or something.
You multiply it by n cubed. Okay.
so, what are the basic ideas? Well this goes back to the 60s.
It's still sort of relevant.
It's less relevant that it was. I'll say a little bit about that.
But, it's still not a bad indicator.
so, basic concept is a FLOP, that's supposed to be a Floating Point Operation.
And roughly, it's something like an addition, a
subtraction, a multiplication, division of 2 floating point numbers.
Now, there's some older definitions where.
It's an addition plus a multiply and you
might include things like square roots in there.
All of that is nonsense because square roots cost like 50 times more.
A division is something like, cost you 20 times more than an add.
So for all of these things, the, the key idea now
is just to not to take the numbers too seriously, right?
They're to be considered to be within a factor of ten.
Right?
And that, that, that's the right number. these,
the flop was invented sort of when, you know,
main CPUs didn't do floating point arithmetic, you know.
Then basically had to like send a message over to
some giant board that was like your floating point board.
And all sorts of stuff happen.
And then when that multiplying was finished.
Right? Something would come back.
Everybody, you know, and so in that setting,
sure, counting the number of flops made total sense.
It's ridiculous now right because, you know it, with the number
of, I mean there's all sorts of things you could do there not even getting into
like minor details but you can do crazy stuff all in parallel and stuff like that.
So, okay.
>> All right.
So, the, the, the idea, and this traces
to the 60s and is still actually somewhat relevant.
it goes something like this.
To estimate the complexity of an
[UNKNOWN]
what you do, is you'd express the number FLOPs.
it usually, it's a polynomial.
It's a polynomial function of the problem dimensions.
And then what you do is you drop. All but the leading terms, right?
So, you know, if it's something like, you know, 1/2 n cubed plus 2 n
squared plus blah, blah, blah - something like that - you just say, well, it's one.
The complexity's 1/3 n cubed.
And honestly, these days, the 1/3 doesn't make any difference at all.
Because you're not shooting for - you know, factor of
2 is ridiculous.
You could not predict using these methods, right?
Factor of 10 is the right, the right amount.
Now you know this is not an
accurate predictor of computation on modern computers.
Right.
It, you can get it within a factor of ten
if, if thing you know in the absence of certain bad
things and actually if things are done really well right
like it's optimized code it once again becomes a reasonable predictor.
So it's kind of an interesting thing.
but, you know it's still, it's, it's very useful as a rough estimate of complexity.
It's still, it's still gives you that, right?
but, you know, do not come along and say oh, you know, this
one is half the number of flops therefore, it must run faster or something.
I mean things like your, your access pattern
and, and, and Locality of reference and memory
are way more important than factors of to
and flop counts, and things like that, so, okay.
I mean another nice thing is that a lot of these things are
done for us by people who know what they're doing, namely those who
write like compilage and things like that, so that, that's the good news.
Okay, so.
let's look at some, operations and this is a conceptual division.
It's a vector-vector, matrix-vector, matrix-matrix product.
And I should say this is something you should probably all know about.
there's no reason not to.
There's something called BLAS and it traces to the 60s, 70s.
It's called it, it stands for Basic Linear Algebra Subroutines.
Right?
And believe it or not, the, the standards of these are written in I, well four tran.
That gives you a rough idea.
OK? And by the way still.
So, and then they have of course, these are
re implimented in all sorts of real languages actually.
So anyway this is basic linear and they call these
things, this is called blass- Level one, this is blast two,
and that's blast three. Right.
You don't need to know that, but actually that's not true.
An educated person just knows these things.
You don't have to know the details of these things.
But you should know but there's blast level
one, blast level two, and blast level three.
And I'll say a little bit about what, what
they mean and what why it's important to know about.
Okay.
So vector-vector operations this would be something
like, you know, calculating an inner product, right?
You calculate an inner product, yeah, you can count the number of flops.
It's like 2n minus 1 flops, 2n if n is large, right?
So that's something like that.
So you'd say it's order 2n or something like that.
Another one would be like a sum.
Something, something like that.
by the way if you're taking an inner product it,
it what, would also depend if, if the vectors are sparse.
Then, needless to say this could be way faster.
And it actually depends, I mean,
quite a bit on You know, the data structure
used to describe sparse vectors and things like that.
So.
Okay.
the other thing let's just say about an inner product, I mean just to give a
hint, sort of, about modern stuff, is something
like an inner product, is, is easily paralyzable.
Right so so on a modern on on some modern platform it'd be
divided up if you do an inter plot of 2 million long vectors
they have multiple threads and each one would be taking inner product in
chunks and then those would be added together in things like that so.
That's.
I'm just saying.
That's by the way how you get the 20 millisecond,
30 millisecond numbers to solve a thousand equations on that laptop.
Okay? So, but it's not bad.
If you want to think of something serially stepping
through like, please get me the next x i.
Please get me the next y i. Please multiply them together.
And then, plus equals that on some, some, you
know, you're inner product plus equals x i y i.
If you want to think of that. Go ahead, it doesn't hurt.
It's charming, it's simple.
But, you should understand that's not how it works.
I mean, unless you write the code right like that, right?
In which, in which case it'll run, whatever, 10 times slower, like so.
Okay. All right.
Next one is matrix vector product.
That's y equals a x.
so here you multiply a matrix by a vector and what is it?
Well, you know there's many ways to think of it.
But, one way if you're just counting FLOPs, it doesn't matter, is that
it's basically, you take the inner product of each row of A with x.
Okay?
So, it's basically calculating M inner products in our N.
And that's something like 2 MM.
Right?
And you forget the 2. This is like,
that's a joke. Right?
It's order M, N. Right.
So, it's the number, it's vaguely number entries in A.
now if A is sparse, this is obviously a whole lot less.
Right? And but actually,
[INAUDIBLE]
how fast it is, is quite a complicated topic and depends on exactly what
data structure's used to represent the sparse
matrix and all sorts of stuff like that.
Right?
But Again, it's not easy at least at a, without looking,
getting into the details to realize that, if a matrix is a
million by a million, with maybe 10 million non-zero entries, you
only have to multiply it by the 10 million non-zero entries, right?
So, okay. and it can be less,
it's a structure. so, for example, here's one.
if I gave you the matrix A.
As a, a low rank product, right? So if I gave you the matrix A.
Let's suppose a is million by million. But I give it to you as a million by ten.
Multiply it by ten by million, right?
So it's what you would get if you ran pca, or sbd, or something
like the low rank approximation.
>> Number 1, you'd actually be able to store it.
You can't store a million by million matrix, right?
So, but you could store it this way.
And the multiply would certainly not be a million squared.
It would be way less.
It's order of a million.
In fact, it's like 20 million times that because, 10 by 10.
Everybody got this?
Right. So, okay any of BLAS level three.
These are matrix matrix products.
So, that means multiply two matrices.
I mean that's the canonical one, is multiply two matrices.
And here, it's actually the product of the dimensions, basically, right?
If these are like m by n, m by p.
It'd be, it's mnp.
So it's the product of the, of the dimensions.
Okay.
that are sparse, it's much less and you get different things
if if they're symmetric, you can divide things by two, I mean,
that typically doesn't matter but that's, that's the idea.
and now, this is probably where I, I mean,
again, we, these are now things We're not going
into these details, but I'll just say a little
bit about it so that you appreciate it, right.
So, everybody here knows how to multiply two matrices, right.
I mean it is, you could write a five-line
C program that would do this, every here could, right.
It can be truly trivial. It could be like how do you get C i k?
Why, you, you walk along A and you walk along
B and multiply the, the associated entries.
And then add it all up.
right?
So, you could write a five line C program. Right?
It'd be a bunch of fours and wouldn't even be five lines, probably.
It'd be 5 lines with a lot of commenting but, okay and you can compile it.
And, of course, it would do, like you know whatever.
Two MNP FLOPs by the time it runs. in fact.
if you were to do, if you were to use
LAPACK, that, that's the, again, this is just for cultural background.
You do not need to know this although it's not bad to know it, okay?
So, LAPACK is just, this is an open source package of linear algebra routines.
and, you know, multiplying two matrices.
This I mean doesn't seem like there be any if you don't think about
this doesn't seem to be any subtlety there right we all know what CIK
is, its sum over J, AIJ be JK right I mean it
don't look like a whole lot of room for optimizing there right.
I can give you something a little
more complicated like sorting or something like that
or finding the shortest path But then you
can imagine somebody smart could do something better.
Everybody, everybody see what I'm saying here?
I mean, come on, calculating the product of two matrices?
You know, one person's implementation, you can beat a naive one by a factor of 100?
That doesn't seem right.
Does anyone know, by the way, how that's done, and why?
Its actually kind of interesting.
I mean, this is just for cultural background, but its good to know.
Here's how you multiply two thousand by thousand matrices.
Here's how you do it if you use this package.
What happens is, thousand by thousand matrices
are broken down into things like ten
by ten block matrices with blocks that are hundred by hundred.
Everybody following this?
Now, matrix multiplication you can do block wise.
Right.
So, what that says is you end up multiplying 200 by 100 matrices
by blocking them and then each block you do 100 by 100 multiply.
Okay.
Now the 100 by, how do you do 100 by 100 multiply?
It's blocked down to 10.
So each of those is a ten by ten matrix of ten by ten matrices.
Everybody got this?
And then so what's really happening is you're
doing giant piles of ten by ten multiplicates, and
by the way the parameters of ten I made up are totally made up and it would differ.
Depends on your architecture. Alright?
Okay, so that's called blocking.
And the idea behind blocking is you can arrange it if you do it just right.
If the block's in the right sizes like that ten by ten
multiplied might hap, might happen in, you know, the ten by ten multiplied
might happen with registers or something, right.
And then the hundred by hundred might fit in L, might just fit in L1, right?
And stuff like that.
So if you do this right, if you do your blocking right.
>> Then, you, and not only that.
You could do part of it in parallel and stuff like that.
If you have multiple cores, you know, somebody's multiplying these.
Somebody's multiplying these.
You exploit all this stuff.
And that's actually how matrices get multiplied.
OK?
On,
at, at least on a, you know, on sort of a modern processor.
Right?
So, that, that, that's the idea.
And, and if you wa, and the memory access pattern is such that there are never
[UNKNOWN].
And by the way, how do they, how do they choose the block sizes?
That's they have something called atlas, right?
That's
[LAUGH]
really kind of cool.
I mean, you don't need to know any of these things, right?
But it's sort of fun.
That's automatically tuned linear algebra software.
this all open source stuff.
And what it does, it takes your computer and it just said, I try 10, right?
And it does it, and it factors matrices for a while
and it's say I'll try blocking it as 9 by 9.
It tries 9, it tries 11, it tries 12, and it just does
bisection or just tries a bunch of them and then half an hour later.
It's decided on your architecture what the
right cash, what the right blocking sizes are.
Okay, completely imperdible, everybody got this?
So, and when you see these curves their absolutely amazing.
It would be things like, you would see
something going like, here's, here's a function of n.
And say it would be, it would be looking like that.
That's a log, log plot.
Then all of a sudden you would see this.
[LAUGH]
Right, and, and you'd say, and by the way, what would this be?
Yeah, so this is something that's just barely fitting in some cache, right?
And this is something that, that.
This is called bad blocking.
This is, this is where your block sizes were exactly You
know a handful of things to big to fit in some cache.
So basically it's a disaster, something like that.
In fact this can even go down that's the cool part.
This can be nonmonotonic but any way okay.
So that's more than you need to know about all of this but it is good to know.
actually it is good to know even
if you're interests are quite theoritical or mathematical.
This is good to know beacuse if you do theorical or very
mathematical things You know, here's, this is how it hits the road.
It hits the road right here, in, in LA Pack in most cases.
Well, except in giant sparse things, but we'll talk about that.
But a lot of times, you do some cool statistics, you do any you do
[UNKNOWN],
it's this.
Now we're going to talk about solving linear
equations, and The one way to do it.
We'll first talk about classes.
classes, of, of, matrices. Where solving Ax equals b is easy.
And then there's a simple trick that puts it all together.
So, one would be this.
diagonal matrices.
Right, so if A is diagonal then solving a equals b.
Is simple.
oh, by the way, one interesting thing is
people will, when people write A inverse b,
mathematically it's what you mean. That's what x is.
but it is universally understood, even if people write it in an
algorithm statement that it doesn't, that you do not parse it this way.
A inverse times b because that would suggest
code like, you know, a.in, you know, times
b and eh, eh, except in extremely rare cases you would never compute it that way.
But,.so, a inverse be is sort of like a, a
phrase that when someone describes an algorithm you should see.
And it doesn't mean form inverse and
multiply by b Mathematically it's the same.
The semantics is identical.
But you should understand that this is so common so, and people will even
say things like inverting that matrix is
easy and they don't mean inverting the matrix.
They mean solving the equation. Okay, so this is just, so, okay.
so here,
you know, computing a inverse b. If, if A is diagonal.
It's incredibly easy. It's just you divide.
It's N FLOPs. You know, roughly.
Whatever it is. OK?
the next one is interesting. Right?
It says that a matrix is lower triangular. Right?
And so, that says, you know it looks like this.
A11, A21, A22, you know, something like that.
Okay?
And you're multiplying this by, you know x1
up to xn and that's supposed to be, you know b1 up to bn.
Okay?
And yeah, how do you, how do you do this? You look at the first row.
The first row says a1, 1, x1 equals b1, so x1 is, is obviously a b1 array 1, 1.
Right?
Oh, and by the way, what exception do you throw if a1, 1 is zero?
What if we have a divide by 0 exception?
Yeah.
You throw a message back which is like A singular.
So, so I'm not going to calculate a A inverse b for you.
Right, so okay.
So once you know x1, right? And there's lots of ways to think of this.
Once you know x1, the second row says A21 x1 plus A22 x2 equals b two right?
But you know x one so you just plug in that value subtract it on
the right and then divide it by A22.
So the algorithim looks osmething like this an
dyou can sort of calculate the number of flops.
I mean it doesn't matter but you can kind of get a rough idea what it is.
>> This thing is kind of growing.
the, oh this, and it's got a, a great name.
This is forward substitution.
And the reason it's called forward substitution is you calculate XI.
And then XI, from then on, you substitute
that value into the equation. Right?
And so, you can see, you do kind of like.
Length 1, length 2, length 3, length 4, and it's growing.
And you know that the sum of those
numbers is something like the square over 2, right?
And if you do it exactly right, it's just like n squared, right?
You count blocks, that doesn't matter. It's n squared blocks.
Okay, and that's called forward substitution.
And there's a, a similar one, if a matrix is upper-triangular.
And what do you imagine it would be, you know, the same thing.
It's like called, it's backwards substitution.
Right?
You start by getting XN first and going backwards.
Okay?
So these are, and these are, just completely simple.
There's nothing complicated here.
this is really interesting when A is sparse.
Because when you have a lower triangular sparse matrix.
I don't know. Go ahead and give me a rough number.
so for a sparse matrix.
1, and a very important attribute is
number of non zeros, sometimes written NNZ, okay?
So I just roughly, if I have, in terms of the number of non zeroes.
How many FLOPS
does it take to do forward or backward
substitution with a lower triangular, upper triangular matrix?
That's sparse, in terms of the number of
[INAUDIBLE]
zeros. What is it?
Yeah, it's winning.
In fact it equals, kind of, to the n and z.
Everybody got that?
Because, basically each of these coefficiants appears once.
Right? And so, it's about linear.
Right?
So, that means if I have a million by million matrice Right?
With let's say 50 million non zeros, that's you know
on average of 50 per non zeros per row or column.
That's 50 million, 50 million flops.
How fast can I do it on there? This rough order magnitude.
You can assume that's between 1 and 10 gig, gig gigaflops something like that.
So the answer is milliseconds.
another one would be orthogonal.
So if a matrix is orthogonal then in fact solving, well
computing a inverse b is the same as well a transpose b.
Right?
So its just a matrix vector multiplied, its ordered n squared in general
but it turns out it actually depends how the orthogonal matrix is stored.
There are other data structures.
Here's a very famous one. And it looks like this.
That's, that's called, that's a reflection.
Right? That's, that's an orthoganol matrix.
And you can see it is an identity minus a rank one matrix.
Right?
We're going to talk a lot about that later today.
but so there that can be done in order n flops.
Right, so very fast. Right.
So, so for example, one, I mean, it's not completely standard, but
it, you would occur often is one data structure for representing an orthogonal
matrix is as a product these things, these reflections here.
And you would actually give them, you would give the u , u1, u2, u3,
u4 and you can see the multiple, multiple actually would be the same to multiply,
multiply by the transpose which is the
inverse, each of those would simply be linear
in that in that data structure, right and if I now want this permutation matrix so
permutation matrix. Well, it's kind of silly.
Right?
If I solve, you know, P x equals b, all I do is permute the interest, right?
And actually by this 1960's version of things.
When you say, how many flops is that, the answer is zero.
There's no floating point in just doing, you're copying.
You're copying I, I guess you, all you're doing is you're copying,
you know, whatever it is eight bytes or, how many bytes is a, a double?
Eight?
So you're copying
[UNKNOWN]
eight bytes from here to there and around like that.
So, and oh, by the way, if that's a hundred
million long vector, and their access pattern is completely random or
whatever This could, you could end up with, this could be
not, it could be 0 flops and actually take some time.
Okay?
[COUGH]
Alright okay, so alright and now
we get to, how do you solve Ax equals b? and, it's a general method.
It's, it's, it's, it's a high level, It's, it's, it's factor solved.
And there are some very important things that actually to know about it.
And some, if you hadn't seen these
stuff before, you'll learn something shortly and it'll
be quite real I mean, it'll be something you should know, and it's not expected.
So the idea is this, there's two
steps, there's factor and there's solve And factor
says you take your coefficent matrix and you
factor it as a product of simple matrices.
By the way, you've seen tons of things like that.
Like SVD, eigon decomposition, eigen vector, you
know, you've seen lots of these things.
Right?
QR factorization.
Alright, so what you do is you factor it as a product of simple matrices.
Sometimes it's two, it just it could be as few
as two and it could be as many as four or five, right.
Actually sometimes it can be a variable number where some of those things
have exotic data structures, like it could be these products of reflections, right.
OK.
So well, one example is you can actually understand the fft this way.
I don't knwo if people know.
If you know what the fft is, great and I'll say what it is.
it's a specific matrix. It's the dfd
matrix, the discreet
[UNKNOWN]
transform.
It's a Vandermon matrix with a bunch of complex numbeers in it.
Do you know what I'm talking about?
Fine. It basically calculates a discreet
[UNKNOWN]
transform.
And there's an interesting thing where you take an order n, matrix like that, and you
can factor it into a product of log
and matrices, which have very specific sparsity patterns.
Okay, and it turns out, if you do this and then
exploit some of the structure, you get an N log N factorization.
Okay, so, alright, let, let's move on.
So, these are things like, you know, diagonal, upper, lower, triangular,
and so on. And then I mean this is totally obvious.
>> What is A inverse B?
Well, I mean, the inverse of a product is just the product of the inverse.
It's in the other order, like this. And then you parse it.
This is extremely important this way, remembering
that you're never computing any of these inverse.
Right? You're, you parse it this way.
Right?
So, you should actually think of that as a method.
A1 inverse that you call on B. so.
That, that's what happened.
So, it, you just un, you just, you unwind, right?
So you write A as a product, calc, calculate
inverse B, you call the inverse methods from backwards, right?
And you just unwind everything. So it looks like, that's the idea, right?
So, so that's the idea. now this immediately has some.
Pretty cool implications. It says the following,
and this, this is like very important.
It says that suppose you need to sol, had solve a bunch
of equations with the same coefficient
matrix but with different right hand sides.
Okay, that comes up in tons of places, right?
I mean Void, tons. Right?
so suppose that's your problem. Right?
And what you, by the way, you'll see tons of cases where this comes up.
Okay. Then, what it says is, you factor
the matrix.
That typically costs more we'll talk about that
in a minute, then the, then the solve.
Right?
so you factor it and you cash the factorization.
Now what that means is it says that if you want, if you have, if you want to
do multiple right hand sides you get one factorization, here, and m solves.
Okay?
So, for the most common methods by the way,
solving linear equations, the factor plus n cubed, the factor steps.
And the back solves cost n squared the solve steps.
Okay, and now I'm going to ask you a question.
To solve 1000 linear equations with 1000 variables on
that machine is, I'll just make up a number.
We'll fix it as, I like your number. 30 miliseconds.
It's 30 milliseconds,
and I want to know how much time does it take to solve two equations.
AX equals B1, AX equals B2, with the same A.
It's 30 milliseconds to do one. How many is it to do two?
By the way, the naive method just calls, solve AB twice.
And the answer is, what did we say, 30 millisecond?
60 milliseconds.
And what is the correct answer? 30.
It's studen, is 30.
OK? All right?
How many
does it. So, by the way, that's pretty weird.
Everyone got that, that the cost of solving two linear equations, if you,
if you, if you cash the factorization, is actually the same as one.
And I don't think this is obvious, frankly.
it is not obvious unless you knew this, right?
It, it's an absolutely real phenomenon and it's, it's the key to tons of things.
Okay? Because of things being sort of either
viable or not. Right?
how much would it cost me to solve, let's say, 25 linear equations.
25 different righthand side same A on that laptop.
I'm doing 25 times more work. And time is the same.
The answer is yes because the factorization dominates everything.
Once you factor you can it simply.
If whatever your background is, you didn't know it, this
can actually have a huge impact on all sorts of stuff
So let's look at some factorizations.
By the way, I'm not going to
cover factorizations here so that's an entirely different
topic and they're not that hard actually
to, to work out how these factorizations go.
It's a whole world in itself.
But this just so you know at the highest level.
We can't go into the details of factorizations.
So. Here's the most famous one.
Is LU factorization and by the
way it's completely equivalent to gaussian elimination.
It says any nonsingular matrix you can factor
as PLU.
Where P is a permutation matrix, L is lower triangular,
U is upper triangular and obviously L and U have
to have non-zeros on their diagonals cause if there was
zero on the diagonal of L or U, what would happen?
Yeah. It's non-invertible.
The determinant is zero and so you have no right solving Ax equals b, okay?
So, okay.
So, and this factorization,
we're not going to go into it although it is
just Gaussian elimination, the cost is order n cubed, FLOPs.
Right?
So, that, that's how that works.
Okay? And so, here it is.
Here's how you solve a set of linear equations by LU factorizations.
So, you're given A and B.
The first thing you do is you factor A as PLU.
That's N cubed FLOPs.
you then, then how do you solve? Well, you go backwards, right?
So,
you solve PZ one equals B, that's like
zero flops or whatever, I mean, whatever you like.
It's, you premuting the entries of B.
Right, then you do forward substitution back substitution, and
you can see this is like two N squared.
Right.
And so, I mean the numbers don't really matter
here but the point is n is typically big enough.
N would have to be tiny, before like 10 or something, then all of this is off anyway.
But n would have
to be, so if n is any reasonable number 10, 100, 1000
or something like that, then the n cubed thing is bigger, right?
And so that's the idea.
And so it's basically about n cubed, right?
It's, it's, it's the, it's, it's the complexity, okay?
So that, that, and that's LU factorization.
the back solve, notice it's order n squared.
So, in fact, you would say the following. the, let's see,
how shall I put it.
Okay, here's a good way to conceptualize it.
one way to think of solving a x equals
b for different right hand sides, is something like this.
The first time you solve it, it costs you n cubed.
Thereafter At, after you ca, you factor A and cash it,
you get a discount by an or, by a factor n.
Everybody got that?
Right?
So, so the first one might take you half a second and then the next
one is going to take you submillisecond, right?
I mean, so these are actually very important things to know about.
Okay.
now, sparse LU factorization. this is quite interesting.
And I'm not going to get into the details,
but this is actually very good stuff to know about.
if any, if you have a sparse matrix A, the most common
factorization method is to factor it into a product of four matrices.
and you know by the way they have other names for these you know p1 and p
2 were sometimes called you know row and column
permutations or something but that is what they do.
That is you factor it into a product of four things a permination LUP okay.
Now the fact is mathematically you
only need one of those permutations to ensure the existence of a factorization.
Right.
That's the mathematically because it follows from the previous one.
Right, the reason you'd use two and reason you would use a permutation both
on the left and on the right is that the L that's actually happening here is this.
As P1 transpose A P2 transpose
Right is L U.
So what you're really doing is permuting the rows and columns
of A, and then calling an unpermuted L U factorization on it.
Okay?
Now, by the way, if you choose P1
and P2 wrong, L U factorization doesn't exist, right?
So let's assume you've chosen them correctly and they work.
It turns out that A is sparse.
The sparsity pattern of l and u
is, depends on the permutations.
Okay so if you pick a good permutation l and u are quite sparse.
And that, of course, is going to affect your solve time.
Because, your what, your solve time. We've already looked at that.
Back and forward substitutions.
Basically n and z. On L, and then Z on U.
Alright.
That's the number of flops.
So the point is, if I pick 1 permutation P1 and P2, so that, L has
a milion non zeroes and I pick another one
so that L has a hundred million non zeroes.
Right.
It's not just a memory issue but also a run time.
Everybody see this?
So okay, this is actually quite fascinating It turns out
that when you do an LU factorization on a matrix, right?
Then the permutation can have a huge effect on the sparsity of the result.
How do you choose P1 and P2?
Well, you can ask lots of questions.
For example, obvious question would be, how
do you choose row and column permutations of
a sparse maker so that it's, it's
LU factorization has the fewest numbers of non-zeros.
It's a totally natural question and it's mp heart.
Okay? So, fine.
so what do people do?
Well, it's, it means everything, these are all heuristics.
Okay?
And, actually, they're more than heuristics, right?
it's, they're even,
they're very stringent requirements, right?
Because what people do who worked on sparse matrices is the following.
If you want to, so what happens is, something like this.
Well, you choose these entries using various methods, right?
One of the methods depends on the sparsity pattern,
it also depends on the values of A, right?
As you're going along you might choose, the P
is actually built up while you're doing the factorization.
Okay?
So that, that that's how this works.
And, and I'm not going to go into any details, but that's sort of the idea.
And it's, it's worth knowing about this.
Okay. And, what's the cost?
Well it's usually, way less than n cubed.
Right.
By the way, if it gets even close to n cubed It means you have
implemented an extremely inefficient method when you should
have just called a,a dense LU pack routine,
right, which would've blocked everything nicely for you.
That's, that's all it means, right?
But, but this is a big deal. so you know, if you're lucky this is
>> This is actually I'd, a, a, a very big
give, by the way, everyone here is a beneficiary of this.
I can tell you that right now.
so, and the reason is very simple. CVX, for example, calls SDPT3 or SEDUMI.
All of these things use sparse matrix methods to actually do the heavy lifting.
That's if you profile any of these solvers, that's all they're doing.
The only thing they're doing is solving linear equations.
So, there's nothing else.
Linear, And they're literally solving a x equals b.
Now a is symmetric, but that's another story, it doesn't matter.
Same idea goes, right.
And they're all spa, and, and they're exploiting sparsity very heavily.
Okay, now they, what's at, there's a little, I mean it's a bit funny, right,
because.
You know, if someone comes along over your shoulder, when you're
plot feeling with CVX, and somebody says, what are you doing.
Well, I'm doing convex optimization, and they go, well, you know, how
do you know it works, and you go, that's the whole point.
It's just like, it's, it's tractable.
I mean, we can do that, I mean, you know,
had a whole course, and that, that was the whole point.
There's other problems you can't solve, but these you can really solve.
The, as a practicle matter the only reason you can solve them, as a practicle matter,
is because the euristics of choosing the permentations
for solving linear systems of equations are so good.
Everybody got that?
Don't, but do not tell anyone outside this class that.
Right? This is a secret we keep together here.
Right?
So, if anyone asks you, what, why does that work?
Is that heuristic?
You say absolutely not.
This is not a heuristic. and in theory it's not,
right.
If, the worse thing that can happen is your
permutation fails spectacularly and you avert to n cubed.
That is not good.
[INAUDIBLE]
10,000 or 50,000. So, but now you know.
This is why, roughly speaking, you can solve equations.
10,000 variables, 100,000 variables, a million variables.
Million by million. You can, kind of, solve, sparse equations.
And, so, the correct way to say it is if the gods, god or gods, who control The
[UNKNOWN]
of permutations used in sparse factorization are smiling on you, and
that depends on the offerings you have made to them likely.
If they are smiling on you then you can solve equations that would just otherwise
be completely and totally out of question, everybody
got this right but there are some sparsity
[UNKNOWN]
Where you don't have to rely on
the gods of heuristic, sparse permutation selection algorithms.
And we'll talk about some of, well I'll tell you one right now.
If the matrix is banded, you don't need to make any offerings at all.
That's, it's really going to work.
Always, m'kay.
But you know for general sparse matrices this is not true, okay, so.
Now this is very
good stuff to know about even if you, if you're more interested
in sort of more theoretical stuff, this is very good stuff to know about.
So now, well we're going to encounter more, I mean what we're going to encounter
are some things that are generally not
just square, but you actually have sine structure.
There's the Cholesky factorization.
If you have a positive definite matrix, A, then you can factor it as it LL
transfer \g, but it's not.
What it says is, it's an LU factorization, with L = U.
Well kind of duh, or aesthetics tells you that right.
because what's a matrix that's symmetric.
Well it's transposes the same, or you swap rows and columns, the same thing.
So, of course, the U and the L, sh, that's all this is, okay.
And L is lower triangular And in
fact, generally people require just to normalize that.
That the entries are positive.
One very interesting fact about Cholesky factorization, is a numerical one.
It actually is stable.
Again we're not really talking about that.
It means you do not have to permute, you, you
can say immediately, that if you calculate that in double precision.
You aren't, you're not going to have any, any trouble, right.
That your er, you can bind the error in L.
And stuff like that.
As a, as a practical matter, it just works.
In this case, and
the cost is 1/3 n cubed. And, you know?
Again, the difference between 2/3 and 1/3
is, like, zippo in, in current things, right?
Where cache, cache misses is way more important than some factor of two.
But, you know?
You get the, rough.
[INAUDIBLE]
, I mean, this, you save because.
Well, you don't have to calculate separately in an L and a U, right?
You're just calculating one L.
Okay.
So, how do you solve things by Cholesky factorization.
by the way, this would come up in least squares, right?
That's your most, that's your most
basic convex optimization problem, least squares.
You can do it by Cholesky factorization.
So.
Here, you'd factor A, as l l transpose,
and you do forward and backward substitution, and
the cost is about 1 3rd n, cubed.
Okay, so that's, that's the, that's the idea, that's Cholesky factorization.
Now sparse Cholesky factorization.
it looks just like sparce LU, but there's actually a
very critical difference, and I'm going to say what it is.
And you probably don't need to know it at
this level, but these are good things to know about.
So here, here it is.
as before, there's P one and P 2, you do pre- and post-multiply.
Or, if you like, you can say that what you, another way to
write this Is, is to say just to think of it this way
first you permute A then call a Cholesky factorization on their perumuted version
right so that, that, that's a a good way to think of it.
Oh and by the way.
Everybody has a different thing, you know.
Sometimes it's lower triangular, sometimes it's upper triangular.
The permutation either goes like
here, or there.
So you always, anytime you actually look, use something that does a factorization.
You have to look at what, exactly what they mean by it.
And it'll probably be different from these notes.
Right, so, but.
That's, this is too bad. There's no standards.
Okay, now here you always have a Cholesky factorization.
But the, the permutation can have a huge effect on the sparsity of L, and therefore
on the memory required to carry this out, and second, and of course, it
also has, since the number of nonzeros in L is basically the flop count.
That's going to affect the run time.
Okay.
Now, unlike the sparse LU factorization where you're actually choosing
P1 and P2 as you go, as you factor the matrix.
Right?
Because you actually looking at entries and
you want to, you don't, want to avoid numbers.
They're too small.
If you wanted to, if you're going to divide by something.
Here it turns out that P can actually be chosen
entirely on the basis of the sparsity pattern of A.
Okay?
So, in other words, all that mat-, you take the sparsity
pattern of A, you don't even have to know what the
entries, the numbers are going to be, Right, if you can make
a graph out of it that's a very common thing, and
then you chose an elimination ordering right but you do this once right.
So, that says if you are going to do something like
solve ax equals b with different values of A multiple times.
When A is sparse definite, you can actually cash the permutation.
I don't know if you got that but this will be,
this will, the usefulness will be clear in about two weeks.
You can cash the factorization.
and that's the idea.
anyway, so the cost is usually you know, way
less than and and all that kind of stuff.
So that's, that's how that works. okay,
another factorization we going to see is the so called LDL transpose factorization
and that's this if you have a symmetric matrix not necessarily not necessarily a
positive definite you can factor it as LDL transpose right and D.
P is permutation matrix.
l is lower triangular unit, so you have l i
equal 0, d is block diagonal, and it's either got
1 by 1 blocks or 2 by 2 to blocks, right?
So and you get sort of the same thing, right?
So this would be kind of one, this would
be one of the methods of, of sort of choice.
Oh, and I think I, I can give you a rough idea of what
>> What backslash doesn't MatLab.
I mean, just roughly. It's something like this.
When you type A backslash B, I think it looks at the diagonal of, of A.
And if they're positive,
or it looks, it quickly looks.
It compares some entries quickly to see if, if
it's possible that A is positive symmetric positive definite.
Right?
So, it guesses that. Then, it attempts the Cholesky factor.
If it finds that it is, if it, if there's some evidence that A is symmetric,
and positive definite, like, for example, it has
a negative on the diagonal, it doesn't do this.
It's all over.
It's not positive definite.
So, it attempts to do a Choloesky
factorization because it's numerically stable and fast.
Okay?
But, if it turns out A wasn't, then the Cholesky fails.
Then, it falls back and does something like an LU or something like that.
I mean this is roughly what it does.
Right? So, okay.
And it won't a sparse factorization unless
A is like explicitly stored flagged as sparse.
Okay. Now we're going to talk about something.
It's about an elimination method, but it actually ties to a whole bunch of stuff.
And it turns out it's actually kind of the
same as selecting the permutation, but we'll get to that.
so what we'll do is let's just take x equals b and block it.
So, all we're going to do is
we're going to, we're going to write it as blocks.
we're going to partition x that'll partition b, and we're going to write it
that way. So, it's just two sets of equations.
It's a1 1x 1 plus a1 2x 2 equals b1, and
then a2 1x 1 plus a2 2x 2 equals b2, right?
So that's, that's it. On this block, you know.
And let's see what happens.
What you do is from the first equation, basically, is a one
one, x one plus a one two, x two equals b one.
And what you do is assumingn a one one is non singular, right?
You put this
over there then premultiply the a one one and
that, that says that you can write x one.
As a function of X2 and you get this, okay.
And so the idea is now we're going
to eliminate X1 from the system of equations, right.
So, all we do, is we plug this value of X1 into this, right.
by the way, this is, this is actually just, sort of, like you know, this,
this substitution stuff.
And with, but with matrices, blocks, right?
So, what happens is, if you plug this equation into that, right.
You get A 2 1 times this thing, plus A 2 2 times x 2 equals b 2.
And when the smoke clears, that's a set of equations for x 2.
And you get.
This equation right here okay.
So this just you've just elim and you
would say things like you have eliminated X1 okay.
This matrix you have seen before that's the shirt compliment
okay and you've seen it in any legitimate class you've taken.
Because a legitimate class uses linear algebra, and then,
you can't use linear algebra and not see that, right?
So you take probability class, that's conditional, that's conditional, you
know, covariance or some, that's
the conditional variance, covariance matrix, right?
You take other things.
Take circuits, mechanics, anything, physics.
This comes up.
Okay, so that's the sure compliment.
so what it says is that when you eleminate variable, a block
or variables, or a block of variables from a set of linear equations.
You're left with a set, with a
linear equation that involves the sure compliment, okay?
So actually, that's, I mean that's good thing to know about, okay.
So, by the way you would form this matrix, you'd solve X2 and then
you would recover X1 using that formula, okay.
And so, here's a completely generic, completely
generic algorithm for solving this block equation.
And it says this It says, yeah I
mean you're going to form, you're going to form this matrix.
that's the sure his thing compliment, it's called S here.
Right, so you, you form this thing, a 1 1 inverse a 1 2.
Now, by the way, that can be formed lots of different Wait oh,
let me ask that. how would you if there were no structure
in a11 and a11 was let's say a 1000 by 1000, and suppose a12 was 1000 by 10.
Please tell me how you would form a 11 inverse, a 12.
What would you do?
So A 1 1. Thousand by thousand.
Inverse times,
A 1 2 is a thousand by 10.
How do you, how do you calculate that product?
It li-, it's a product.
How do you do it? A 1 1 cash.
And do 10 back solves. Done.
Perfect. Okay?
So. All right.
So, but, you just write, people just write stuff like form A11 inverse A12.
Okay.
And assuming you know what you're doing. Okay.
so, you form this, your compliment.
then you solve this, your compliment equation, and
then, you use this equation to reconstruct X1.
Okay?
So, that's the method. All right?
And, you know let's, let's, let's figure out how much it costs.
And the answer is, it's F plus N2S.
That's the cost of factoring A11.
And S is the solve step. Right?
these are the back solves. Right?
As you, as you said, that there's one factorization and some solves.
step 2 cost this.
And what we're going to do here is we're just going to
assume on the second step as dense.
I mean, because you could work it out, when it's not Right and then the
the the last step is a something that's again we're going to just assume its
dense generic nothing and we'll just do
we'll just solve that equation and that's going to
be a N, N, N two cubed right so the total thing looks something like this.
Now.
Here's what happens. If, if you apply this method now, to a
matrix that has no structure whatsoever, then f equals this, and s equals that.
And you plug it into this formula, and you get back, actually, exactly.
The fact you get, there's no savings what so ever okay?
So that was a completely silly thing to do.
I mean, it doesn't matter that's how you do it right.
Actually, it's interesting because it shows you,
it shows you what blocking is right?
because there is, this is blocking 2 by 2 now you can imagine what blocking
10 by 10 looks like, it's just this. Okay, so that's it.
Here is where it's useful.
>> When A11, has structure, right? So, the factor, with the factor cost and
the solve cost on A11 is low, then, this delivers, this delivers big savings.
Here's an example. Suppose that A11 was diagonal.
Right? That means the factor cost
is zero and the solve cost is n1. That's just, just divide.
Right?
Then the number of flops is this, and this is absolutely amazing,
so this is that if you see a set of linear equations
that look like this. Right?
So, It looks like that.
I'm going to, I'm going to draw the, non-zeroes.
Right?
Oh sometimes people call this a spy diagram,
I have absolutely no idea why, but anyway, whatever.
I, it doesn't matter.
I'm just showing you the non-zeroes.
Okay?
Now there's many other ways to understand this,
and and I'll, I'll show A, a couple and
the once the way they are purely
mathematical you probably have seen before, OK.
So and also introduces a very interesting idea so you what that is.
So lets , lets takes this one suppose you wanted to solve A plus BC and
the idea here is that B is like you know it, its so this is low rank.
Looks like that.
And here's, and here's a. and let's assume a has structure.
So let's make a banded. You know?
That's, 'cuz it's pretty, and it has structure.
So, there's your, your sparsity pattern.
That's a plus bc, okay?
it could be anything else. It could be sparse.
It could be.
Forty eight, you know, DFT, it could be anything, you know,
it could be something that you have a fast transformer for.
It could be some wave lit thing or something, I don't care, right?
And then so the people would call this is like structured plus low rank.
That's what this is.
So it could be sparse plus low rank matrix or something.
Okay.
Now the dif, the sparsity pattern of A plus BC is generally full.
Severy, you know, cause, you take an outer,
you just take B and C to be columns.
Or, you know, B is a column and C is a row vector.
Then it's just an outer product, it sprays entrys all over the place.
Right, so, so then, a generic, a generic method there will be no sparsity, nothing.
So, how you can solve these is very clever, and
in fact the first step, there's a name for it.
It's actually called unelimination.
So elimination means you take your variables
and you look at one block and you
say, I'm going to express that in terms of the others and then get rid of it.
Unelimination means that you're going to introduce a new, a new one.
So what you do is you look at (A + BC)x =
b, you write that as A + By. >> Right?
Equals B. Y is actually CX.
Everybody got that?
That's just, of course, they're the same thing.
They're identical.
But, you rewrite these equations like that.
Right?
So, there's, there's, and , and, here, so, you've, i mean you've done unelimiation.
You've introduced new variables and new equations.
Everybody got this?
And you see this.
Now, you're not trained yet but you will be.
And it, once you're trained, what your eye,
what you're searching for is the
following: easily inverted matrices on the diagonal.
And your eye should, once you're trained, go immediately to that minus I.
You will go, it will go immediately there.
Because in the world of easily inverted matrices,
I is up there and -I is up there.
Okay? Everybody cool?
Alright,
now.
I need, see, I just have this fantastic, I just learned about block elimination.
There's a, I'm going to elimi- I eliminate that block.
Now if you eliminate that block you're eliminating y and where do you end up?
Right back here.
Everybody got it?
So, well duh, I mean you un eliminated and
then eliminated the thing you had just introduced right?
OK.
Fine. So the key, is to, is to control.
>> Your urge.
Primal matrix structure factorization urge, to eliminate that minus I.
And instead, you take the system and you eliminate A.
And, in fact, that's not a bad idea because
the assumption here is that A is easy to invert.
Right?
So, A is 48 transform or some wave lifting.
I don't know what it is. It doesn't matter what it is.
Right?
So, that's what it is. OK.
If you eliminate A, so, if you uneliminate.
If you go from this. To this, and then you eliminate a.
Not the minus i.
You get this equation.
And you will recognize that as the Schur complement here, right?
It's like, you now, this thing minus that times the inverse of that times that.
OK. Or other way around.
Some, that's that's Schur complement.
OK? And you get this.
Now, then you solve for ax, and you can, once you get y,
you get x, okay?
And you can actually write this out, just
purely mathematically, and it comes out to be this.
A plus b, c inverse, is a inverse minus, and then this big old thing over here.
Okay, it's a very thing, it's got all sorts of names.
matrix inversion lemma I think is a nice one, and you'll just
do a quick example and then, and then we'll quit for today.
so if A were diagonal and B and C were debts then
that's a beautiful example, right?
So that's called a diagonal plus low rank Matrix, right?
Perfect, perfect example of this structure, right?
then it says you can actually solve it in here
if you, you do the dumb method it's like n cubed.
If you'd used this method which is an un-elimination and an elimination it
actually ends up it goes from cubic in n to linear in n.
OK?
So my claim is, you just, these are not
small differences, right, these are not theoretical things, right?
The difference between knowing this and not knowing it means you can look
at something, and it could be a million by million set of equations.
But if it's diagonal plus low right, now you know.
I can solve that in probably in like sub-millisecond.
But you have to know what you're doing.