Optimal computation
CC Attribution 3.0 Germany:
You are free to use, adapt and copy, distribute and transmit the work or content in adapted or unchanged form for any legal purpose as long as the work is attributed to the author in the manner specified by the author or licensor.
DeVore, Ronald
A large portion of computation is concerned with approximating a function u. Typically, there are many ways to proceed with such an approximation leading to a variety of algorithms. We address the question of how we should evaluate such algorithms and compare them. In particular, when can we say that a particular algorithm is optimal or near optimal? We shall base our analysis on the approximation error that is achieved with a given (computational or information) budget n. We shall see that the formulation of optimal algorithms depends to a large extent on the context of the problem. For example, numerically approximating the solution to a PDE is different from approximating a signal or image (for the purposes of compression).
 It gives me great pleasure and it's an honor to introduce Ron DeVore, who over the last few decades has been pioneering and introducing a variety of methods of functional approximation and approximation theory in all of its ramification. Ron's contribution covers a broad range of issues and problems, deep analytical issues related to methodology for functional approximation. You will hear sort of the outcome and import of this kind of work in his talk later, but let me just say that in our day functional approximation is sort of the key to the digital age in some sense. Machine learning and functional regression, it's sort of the main straightforward elementary question that everybody asks, but nobody knows how to answer, and hopefully we'll learn something in the next few minutes. Later aside, I just learned a few minutes ago that we are actually mathematically related. He's my nephew, mathematical nephew. We graduated under Bojanic, who is my mathematical brother. We are both students. We were both students of Karamata. So the world is small in mathematics and beautiful. So thank you, Uncle Ralfi. Ralfi. I think we're all aware that most scientific problems can't be solved exactly, so we have to resort to some sort of approximation or numerical computation. So we create algorithms to try to resolve these problems, and my talk is concerned with is there a way to tell whether an algorithm is optimal or near optimal? This would be not only a wonderful result mathematically, but would lead us in the right directions of how to proceed computationally. So to answer that question, of course we're going to have to describe exactly what we mean by optimal. But I want to put another caveat into this. You may think of an algorithm that sounds very good to you in theory, but you have to remember that whenever we implement such an algorithm, we do it on a computer, and we're subject to machine errors and even errors in entering the data. So it's not just a question of designing an algorithm, but the algorithm should be numerically stable, and we'll see that ingredient coming up in my talk. So there are many settings that I could describe this, and in fact with my colleagues, especially with Albert Cohen and Wolfgang Dahman, we've studied this problem of optimal algorithms in many different settings. But after some thought about what I should do here in this talk, I've decided to speak about the, in the setting, this clicker is not working very well, so I think I'll stand here maybe to do better. I decided to concentrate on an area which is called compressed sensing. This is a relatively new area in mathematics. It's getting a lot of interest these days. And the motivation for this is the problem of sensing real world signals. So the usual paradigm for signals is that they're band limited functions, and the way to sample them is through what's called Shannon sampling theorem or Nyquist sampling. The unfortunate state of affairs is that many of the real world signals that we're interested in are broad banded. This means that their Fourier transform has large support, and that this type of sampling cannot be implemented in practice. So compressed sensing is looking for a way around this to find another way, an innovative way of sampling that defeats the roadblock of Shannon sampling, and samples signals more at the information content in the signal rather than at the, determined by its bandwidth. So my reason for choosing this subject is there are three reasons. One is that it's very audience friendly as you'll see. It'll be easy for you to grasp the main ideas that are here because it's basically going to be an exercise in finite dimensional geometry and linear algebra. The second reason is that the subject itself interfaces a lot of areas of mathematics, most notably probability theory, theoretical computer science, and functional analysis. But the real reason is that by discussing this subject, I'll be able to introduce the types of optimality that are important in numerical computation and all the nuances. As we'll see, I will introduce three types of optimality when we discuss this problem, and it's sort of the full range that I know of that occurs in numerical problems. Okay, so let's begin. Well, really in compressed sensing, we're interested in analog signals, but the subject is fully developed only for discrete signals. So I'm going to stick to discrete signals until the very end when I make a few comments. So our problem is the following. We have a vector x and rn, and here capital N is large. So for example, if we were dealing with images, n would be the number of pixels, it would be of the size one million, or even more. And the problem is that the game we're going to play is we're allowed to ask small n questions about this signal. And these should be non-adaptive, by which I mean I don't want you to have that the one hundredth question depended on the answers to the previous 99. You must set the questions down in advance. I'm even more restrictive because what I mean by a question is going to be very specific. A question to me is that it's going to be the inner product of this vector x with some other vector. So you can choose the vector v, and that will give you then by taking the inner product of v with x, you'll get the answer to your question. So question is very specific in this case. It's linear information. Now, we can always represent this kind of a system by an n by n matrix phi. The rows of this matrix would be these vectors v that we're taking the inner product with respect to. And capital N is the length of the vector x. So typically, capital N is very large, and small n is very small. So we have short and fat matrices phi. So you see in this world, the Brazilian diets, this is not applied. We're interested in short, fat things. Okay. My question is, can we say what are the best matrices to use in such a system? Or if not best, near best, or good matrices? Now here, good is going to mean what? It means that once we extracted this information y, which is phi of x, so x appeared before us. We asked these questions. We get this vector y, and now x disappears from our site, and all we can use is y. And what we're interested in is how well can we recover x from this information y? How well can we approximate x? This is roughly speaking what we mean by good. Now there's going to be two aspects to this problem. One is, does y contain enough information about x? And the second is, how do we extract this information from y? That part is called the decoding, and both are serious issues. Okay so let's think about what the problem really is here. So we have matrix phi that's mapping a large dimensional space, Rn, into a small dimensional space. So naturally, a lot of vectors are mapped into the same image. In particular, if you look at all the vectors mapped into zero, the null space, this is a large space. It will have dimension capital N minus little n. Since capital N is big and little n is small, this is a large dimensional space. We have a lot of collapsing going on. If you took any vector y and asked which vectors are mapped into y, they can be described as just a hyperplane where you take any one of the vectors mapped into y and then add the null space to it. This gives you a hyperplane, and this entire hyperplane has the same information y. So these hyperplanes which are denoted here in yellow, so this is a given hyperplane for a given y, all the points on this hyperplane are mapped into the same information vector y. Right, so there's tremendous collapsing here going on. Now a decoder, what is it going to do? It's going to take y and it's going to map us back into the big space, rn. So after it does that, we get our approximation x bar to our original x. Now at this point, you should be very pessimistic that we're going to do something here because of this collapsing. So this vector x bar would have to serve as an approximation for all the vectors in this big hyperplane that was mapped into y because none of those are distinguished in any way by this encoding decoding. They're all treated the same way. Well indeed if we had to work in complete generality, we wouldn't be able to say much about this. However, our understanding of the world of real signals is that they have some structure and we want to utilize this structure or take advantage of this structure in going forward. So our paradigm or the way we think of real world signals is that they're sparse in some sense, that there's a set of building blocks and if you represent this signal with respect to this set of building blocks, it either has a very small representation or it can be approximated well by a few terms. Now building blocks, in my case, I'm going to restrict to the case where the building blocks are a basis for our end. You could think of this in other settings as well. In some problems, we may know the basis in which the vector is sparse. In other problems, we might not know the basis. So this is an important consideration. I'm going to actually proceed for the moment assuming that the basis is known to us and then later I'll say something about what we do when the basis is not known. So here's going to be my first way of measuring optimality. So I'm interested in telling you that when we tackle a problem like this, let's try to make a precise mathematical problem to see what we mean by a best or near best algorithm. So here's my first notion of optimality. I want to deal, I'm going to assume that my signal x is sparse. And I'm going to, for now, I can assume that the basis in which it sparses a canonical basis. This is not, I'm not pushing anything under the rug because you can always make a basis transform to get to this case. So what do I mean by a sparse vector? I mean that it has few non-zero components. By let's, given a vector x, the support is the number of i's, the indices for which the vector is non-zero. And I form this class, sigma k, which are all vectors whose support is less than or equal to k. So this is a non-linear set. I mean it's not a linear space. It's a union of hyperplanes. You can take, for example, take any index at t of size k and look at the set of vectors supported on t. And then you union these. So you can visualize this as you have a lot of hyperplanes and the union of these is the set sigma k. So I'm going to assume my signal is in sigma k. And I'm going to ask, given that information, how can you best design an encoding-decoding scheme? So here's the problem. No one is fixed, little k is fixed. And I ask, what's the smallest n such that I can build a matrix phi and a decoder that will capture every one of those guys exactly? And I'm going to tell you right off hand what the answer is. The answer is n equals 2k. Now this is sort of morally justified if you think about it. Why? Because a vector in sigma k is determined by 2k pieces of information. You need to know the k positions where it's non-zero. And then you need to know the value of x at these k positions. That's 2k pieces of information. So one would think you can do this with 2k samples. But of course, my definition of question was very specific. I have to take inner products. So it's not absolutely clear that I can do this. But as we will see, we can. So to describe what matrices, so my goal now is to describe to you what are the good matrices for this problem? What matrices will work? So I want to think of my matrix as built from its columns. So this matrix has columns v1 through v capital N. And I'm going to introduce a property of the matrix that if the matrix has this property, it will allow a solution to this problem. And to describe this property, what I look at is sub matrices of phi. And the sub matrices I get is I take any t columns, any set t of indices, column indices, and I look at the sub matrix gotten from phi by just keeping those columns and throwing away everything else. I'll denote this by phi t. And I will use this matrix frequently. And associated to that matrix is this Gramian matrix, phi t transpose phi t, which is the matrix of inner products. This matrix is a square matrix. It's symmetric and it's non-negative definite. And now here's the theorem. So suppose you have an N. And you want to check whether a matrix of this size will capture every k sparse vector exactly. I claim this will be the case if and only if one of these four conditions hold and they're all equivalent. So the first condition is that the intersection of sigma k of these three conditions, since this is the requirement, that we reproduce all signals exactly. So the first condition is that sigma 2k intersects in all spaces zero. Now notice I'm trying to recover k sparse signals, but here the condition is 2k. Now why do we have that? It's very intuitive. We need to avoid the case where two signals of length k are both mapped to the same y, because then we couldn't find which one it was. While two signals are mapped to the same y, I can take their difference and that'll be mapped to zero. Well that's something in an all space, but their difference will typically be a signal of length 2k. So that's how I get to this condition. And equivalent conditions are the following and these are the simple linear algebra. Another way to state this is to say that these matrices phi t that I form, no matter what selection t of columns I pick, as long as I pick 2k columns, that this matrix should have rank 2k. Another way to say it is that this Gramian matrix should be invertible. So if I can find a matrix phi with these properties, then I know it will recover two k sparse signals exactly. So the question becomes, can I construct such matrices? Can I write some down? Maybe I'm stating something you can't find such matrices. Well what do we need to do? We need to create vectors in r2k. I'm claiming I can do this for n equals 2k, right? But I need to be able to create a lot of vectors because capital N is arbitrary. It can be arbitrarily large. I need to create a lot of vectors in r2k and I need that whenever I take 2k of them, they're more or less independent. They have to be linearly independent. Well here's a quick and easy solution to this. Take the Van der Mann determinant. Take distinct points x1 through xn and form the matrix where in each column you have the successive powers of the given xj's. Then we know that if we pick a sub-matrix consisting of 2k columns of this, we get a square matrix which corresponds to polynomial interpolation. This matrix is invertible. So this matrix satisfies the properties I wanted. So it solves the problem. But I haven't discussed how I'm going to do decoding. Now you can think of many ways to do decoding. I'm going to just put down a naive decoder. A naive decoder would be I'm looking, now I'd like my decoder to be defined for every y, not just for the y's that are images of k-sparks vectors because in practice I won't have that information. So I want my decoder to be defined for every y. So here's something you could do. Given y, you could look over this set sigma k and try to find a z that best fits y, that is that y minus phi of z is small as possible. This is a least squares problem of trying to find a least squares fit to y from the image of sigma k under your map phi. You can decouple this problem. You may look at these spaces that we have introduced before where t is a set of column indices of size k and we solve the least squares problem on this set t. And of course we know we can actually write down what the solution is. This is more Penrose inverse where we have to apply this Gramian matrix, the inverse of it, to phi ty. So this already shows you where we need somehow that this matrix is invertible. We need non-singularity. And then we just search over all these t and find the t star that gives us the smallest residual. In the case that our vector was the image of a k-sparks vector, there will be a t for which this is zero. And it will be a unique t and we'll pick up that t by this decoding process. Okay, sounds great. We solved our first problem. Can we now go home and have a glass of cognac and feel happy that we've done something positive? Well, I'm sorry to tell you that there's bad news here. None of us will be alive when this decoding is finished. Okay? So this decoding that I told you, this naive decoding is horrible. If you had a problem where n is like a million, like an image processing and little n is like a thousand, you can imagine you're going to have to take capital n things taken little n at a time to check. But a more serious problem, it turns out, is that the decoding is not stable. Now you can actually fix the first problem. You could find another matrix. For example, you could take the discrete Fourier matrix, first n rows of a discrete Fourier matrix, and then by some clever decoding, you can do this whole problem in this number of operations, which is a reasonable number of operations for us. So this isn't the real problem that the decoding takes forever. But the fact that the decoding is unstable is a real problem. And let me just explain this to you that we're really stuck here in some sense, because suppose I'm a generous man. Suppose I tell you where the vector x, what its support is. Suppose I tell you, OK, it's supported on t. Now you only have to find, you don't have to find these indices of support. You only have to find the vector x. How would you do this? Well, you have to apply this inverse operator. And this, of course, we know in the examples that I created that this is very poorly behaved. And in fact, now, if I stick to wanting to have only 2k samples, I'm stuck. I cannot do this in a stable way. OK, so what seemed like an easy mathematical solution, a lesson to us, is what seems mathematically OK may not be numerically reasonable. Well this brings us to actually where compress sensing began. I'm not telling you this in chronological order. Compress sensing began about three or four years ago in the work of Immanuel Kandez, Terry Tao, sometimes Justin Romberg, and in a parallel development by Dave Donohal. And what they do is they show how to create matrices that not only will recover 2k sparse vectors, but do it stably. To do this, one has to pay a little bit. One has to increase the number of questions, and we'll see how that plays out. So they make, and I'm going to talk about the Kandez-Tao development, which I think is very elegant and very much to the point, they make two important contributions. First, they're going to tell us what are the good matrices, or they're going to tell us a sufficient condition that a matrix be good. And secondly, they're going to handle this problem of decoding. They're going to show us how we can do fast decoding. So how does this play out? So the property that they're going to introduce for the matrices is the following. Remember these matrices, phi t and these Gramian matrices, they're going to require, we say the matrix satisfies a restricted isometry property of order k. If whenever we pick k column indices, the eigenvalues of this matrix are in 1 minus delta 1 plus delta, and here delta is between 0 and 1. So the important thing is the eigenvalues are staying away from 0, and they're staying bounded. So this Gramian matrix now is not only invertible, but it's nicely conditioned. So that's the new ingredient. You could restate this property as follows. You could just say that what you need is that whenever you apply phi to any vector in sigma k, the resulting vector has norm roughly the same as the norm of x. That's where the word restricted isometry comes in. If we could take delta in this equal to 0 in this formulation, then it would mean that we have an isometry. But we give up a little bit, and we have to give up a little bit. Now what about the decoding? They suggest for decoding that once use L1 minimization in this way. So suppose we have our vector y. To decode what I'm going to do is look at all x which map into y, and then I'm going to pick the one which has minimum L1 norm. This has a very simple geometrical description. If I think of y, it determines this hyperplane, and this is the ball in little L1. I take this ball very small about the origin, and I blow it up like a balloon until it hits this hyperplane. When it hits this hyperplane, that would be the solution. And generically, it usually happens that it hits at a unique point. Now here's the main result. By the way, an important part of L1 minimization is it can be solved by linear programming. So we have a modest investment in the decoding. For example, we always know we can decode with order capital and cubed operations. And here's the main result of Candidaz and Tal. They say if you have a matrix phi that satisfies this restricted isometry property, now for 3k, before I had 2k, probably you could replace 3k by 2k and everything if one studied this enough. But once you have this property, then you can find a matrix with this. If you have a matrix with this property, then you use L1 decoding, you'll capture every x exactly. That's in sigma k. Every sparse guy will be captured exactly, but moreover, the decoding is stable. So now our problem is, OK, we have a new matrix problem. Then we find matrices phi that satisfy this restricted isometry property. And how do we do that? So how do we build good matrices? Well, I claim that for given n and n, we construct such matrices. So we'll satisfy restricted isometry property of order k. For every k less than n divided by log n. Now all of a sudden you see a log n appear in this, which hadn't appeared before. This is the price we pay for stability. You can actually show that you can't discard this log n. It must be there. If you want this restricted isometry property, you're stuck with this log n. So we give up a little. Before we could take k equal to n over 2. Now we can't quite take k at large, but only sacrifice by this logarithm. Now, the question I want to discuss next in the next few minutes is, how can we construct such a phi? What do we want to do? Well, we're back to our problem of creating a matrix with a lot of columns, matrix and our vectors and our n, and a lot of them. Now we want them not only to just be that whenever we pick 2k of them or k of them that they're linearly independent, we want them to be very far from being dependent. We want them to be almost orthogonal in some sense. I'm going to mention to you now three ways you can do this. All of these constructions are probabilistic. The first way is we could look at the unit sphere and our n and to choose our columns we simply pick at random according to the uniform distribution. Independently at random we pick our columns. This would work. Another way is we could create a stochastic matrix. What we could do is take one of your favorite distributions, let's say the Gaussian distribution with mean zero and variance one over n. And for the first entry just make a draw of this distribution and whatever number you get stick it in as your first entry in the matrix. And then independently repeat this every time you want an entry and thereby fill out your matrix. You could do the same thing with a Bernoulli distribution with plus or minus ones or zero and ones. Only one of these constructions will work and will give you a matrix phi with high probability it will work. That is you do this realization then you can prove that with overwhelming probability the resulting matrix will satisfy the RIP property for this range that I advertise and I tell you this is the best range you possibly can have. So we have a way to construct these matrices. I want to make one point very clear here where probability is coming in because sometimes people get confused on this point. We use probability to construct this matrix phi to prove the existence of phi. What we know is we do this construction and almost every matrix we get is going to have the property we want. Once we have a matrix in our hands the algorithm is constructed. There's nothing probabilistic in the algorithm. We apply phi, we apply L1 minimization and so on. Probability is only used to construct this. Because this RIP property is so important, it's of interest to have nice methods to verify RIP. Typically one would think that this is an eigenvalue problem and I'll approach it through eigenvalues. What I want to point out here is that actually there's an elementary track to verify the RIP property in the settings that I just described, for example in the Gaussian or Bernoulli case. So suppose I generate a stochastic matrix by taking draws of some random variable R and fill out the matrix in this way. Then the first thing to observe is you can just write out what phi omega of X is and expand it and take the norm and you'll see that the expectation of this vector is norm of X squared. So I'm going to try to draw out what distributions, what probability distributions would work in creating these matrices that satisfy the RIP property. So this is always going to hold for any distribution, but we're going to need more. What we need to have more is that when we look at a given X and we look at phi omega of X that it's with high probability close to this mean. So we need a concentration of measure inequality, namely that for any given X we have this estimate on the probability that these two things differ by a delta times the norm of X. This concentration of measure inequality will be satisfied for Gaussian, Bernoulli, many distributions and it's usually an easy matter to show that this property is satisfied. And I claim that whenever you have this, then the resulting matrix that you create from this random process will satisfy the RIP property for the range of K that you want. And maybe it's not so important that this is true, but what I want to emphasize is the proof because there are other ways, as I mentioned through Eigenvalue analysis, to prove this result, but there's actually this has a very simple analytic proof that I want to just quickly draw upon. So how would I prove this? What I do is, see I need to verify on sigma K, I need to verify this restricted isometry property which meant that the norm of phi applied to X is like the norm of X for all X and sigma K. That's what I need. So what I do is first I cover the unit sphere in sigma K to some tolerance, let's say delta over 4. Then I know that with this covering, this net of points, for each point in the net I know this, this is a concentration of measure inequality, told me this. With high probability this is true. And since it's with high probability I can apply a union bound and get that as uniformly true for all Q and P. And of course you've got to be careful of how many points do you throw in there, but do this. And then finally, you need to extend it to every vector in sigma K, not just the points in this cover P. And to do that is quite trivial. For example, the upper estimate, you can estimate norm of phi by saying, well I'll approximate X by a Q from my net. And then I have this trivial estimate for norm of phi of X. And the second one is under control because of this. And the first one, well any matrix on a finite dimensional space is bounded, so there is some bound. But of course I won't get the best result this way with M. But then I can bootstrap this. Because once I get a bound for M I can repeat the process. And if you do this you'll arrive at this theorem. Okay, now I talked in the beginning about, you know, the basis is known to me. What can we do if the basis is not known? And people in this subject talk frequently about this concept of universality. Well, these random things, they work for any basis. But I caution you to be a little bit careful on this point. Let's make this precise. To make this precise, if I start with any collection of bases, and if the number of bases in this collection is controlled by e to the cn, n is this, you know, remember the number of rows in a matrix. Then this same argument I gave probability argument can be used to show that you can prove the existence of one matrix phi that will satisfy RIP with respect to every one of these bases simultaneously for the hope, the range of k that you would like, the biggest range of k. The problem, so this is a sense that you can say universality. You can say that if you draw at random a matrix with high, and if you have these bases known in advance to you, this collection of bases, and with high probability, this will uniformly handle and encode well, spar signals with respect to any one of these bases. The problem is that we would need to know the basis to decode. This is good for encoding, but not good for decoding. Well, some people talk about, well, I'll encode and then 20 years later I'll decode. You know, I'll put it in a time capsule. I don't know too many applications in which people want to wait 20 years, but. Okay, that was my first example of optimality, and it was the simplest. Now I want to go to more general. So we know that real world signals aren't typically going to be spars in the sense that I mentioned to you. They're not going to have support that's only a finite number of terms or a small number of terms. They're typically going to have a lot of entries in them. But in general, most of these entries are going to be small. And our idea, remember as we began, was that we said that we view signals as being approximated well by spar signals. That is, they have only a few components that are essential to capture them, and we only need to find those components. So to make this meaningful discussion, I want to talk about how a signal X can be approximated by spar signals. So I introduce this measure here of error. So you can choose whatever way you want to measure error. The norm X maybe for this talk just take X is the least squares L2 and and then given a vector X, we look at how well it can be approximated by vectors which have support K. So a vector is nice if there's a vector of length 100, support 100 that approximates it well. And we think real world signals are nice. So mathematically, we're interested in classes of signals in which we know that they can be approximated well. And typical way to get such classes are the following. You could take the unit ball of an LPN. So for example, if you want to measure error in LQN, and I suggest you to take Q equals 2 here, then any vector in LPN will be approximated by K terms to this accuracy. Now look at this a minute. As P gets small, now here small P is better than large P. Small P means that the vector X, the components of it die off very fast if you order them in size. And you'll have this decay, K to the 1 over Q minus 1 over P. So this can be very, very fast if P is very small. For example, if Q is 2 and you take a vector in a unit ball of L1, then it can be approximated by K terms to accuracy 1 over square root of K. This is easy to see and I don't stop to prove that for you. So what would I mean by optimality in the second sense? Optimality now will mean I have a class of signals, capital K, and I have some norm in which I'm going to measure distortion or error. And I say that the sensing system is optimal if it performs like this min-max. So this min-max says, OK, let's look at how well we can possibly perform. This is the best performance we can make by systems of size n. That is, the matrix has a dimension, the little n by capital N. We take such a matrix and a decoder and we apply this to a given X and we achieve a certain error. And now we look at the super, the errors over all X in the class K. So this is a usual min-max way of doing things. We look at the worst error and then we look at the system, the phi and delta, that give us the best worst error. So we take the infimum over all these. And now I would say a system is near optimal if I have this kind of an inequality with some absolute constant C0. If I had C0 equal 1, I would have optimality. There would be nothing better, but usually you have to give up a little something. So this is interesting to understand for a given class of signals, how does this en behave? What's the target? What's the best we could possibly do? That's the usefulness of en. It tells us the best we could do. And the story is that this, the behavior of enk for classical sets K in our n is known. It was not only known, but it was done a long time ago in the 1970s and early 80s in the context of finite dimensional geometry or geometry of Banach spaces. The main players were Gaussian and Gluck skin. And I just mentioned one result. For example, if I look at the unit ball in L1n and ask what's the best encoding decoding I can do with that in L2 when I measure the error in L2n, then here's the estimates. Notice the estimate from above and below are the same order. The only thing that differs are these constants, C0 and C1. The estimate is basically 1 over square root of n. 1 over square root of n would be what we could do if we chose the n best terms of the vector. And we pay a little price here, this log of N over little n again. Remember we had that before. And this price we know is necessary to pay it because we have a lower estimate. Now I'd like to say so why was this not a solution to the compressed sensing problem back in the 70s? The thing that these players did not look at was practical decoding or how can you implement this encoder? In fact, they never touched the problem of decoding. And a contribution of Candidaz and Tao and Donohull is that they show that if you have a matrix with the restricted isometry property and you add L1 minimization, then it will give near optimal performance. It'll achieve this type of performance. Now not for every Q and P, but for certain Q and P. For example, if you want to measure error in L2, then for every P less than or equal to 1 it'll be optimal. So they give a practical scheme and once you have this scheme, now you're happy because you know that you're doing the best you possibly can. Let me say just a couple words about this because there's, you know, I wouldn't say a controversy but I talk to people in functional analysis and they say, look, these guys in compressed sensing, they're just doing our work over again and so on. This is not a fair statement because as I mentioned to you, they did not discuss the practicalities of the decoding. So their results are useful and interesting but didn't discuss practical decoding. What is the most useful about the results is the lower bound, it turns out, because it tells us we can't do any better. But it doesn't give us real practical schemes to do the encoding and decoding. So I think both camps have made a great contribution. Now I want to go to my third and last way of measuring optimality. And it, to me, if you were in the States, you would call this the Cadillac of all ways of measuring optimality. And Germany would be, I guess, the Mercedes Benz. I don't know what it is in Spain. But anyway, what I want to have optimality in a very strong sense, here's the sense. I'll say that this encoding-decoding scheme is instance-optimal of order K. If the following holds it, every vector I choose. I don't need it in any class or anything like that. Take any vector. When I encode and decode, it performs like K term approximation. Well, this would be terrific if you can do this for large values of K. For example, if you have this and you can return to your first problem, and if you take a vector which is K sparse, this will be zero. So you'll encode, you'll get it exactly. It will also include all the results on optimality for classes. So this to me is like, you know, if you could get this, it's like a dream, but if you could get this, this would be the best possible result. And the problem is, let's suppose you fix the size of the problem, the size of the system you're going to use, the number of rows you're going to use in the matrix. Then the question is, how big can you take K and achieve this instance-optimality? And with L Bear Cohen and Wolfgang Dahman, we've solved this problem for all, if you want to measure the error in any LQ space as long as Q is between one and two. And I'm just going to mention a couple of these results. So there are two things to say. One, the good news. The good news is that you have instance-optimality in L1, and you have it for K of the range that you're the most you could expect, namely K less and log n. Moreover you can build this system and use L1 minimization as a decoder so you can actually build a stable system and have this instance-optimality. And this result, although we prove it in our paper, you can, it essentially follows from the work of Candez and Tao. Now that was the good news. Now here's the bad news. The bad news is that if you ask for instance-optimality in L2, you want to measure error in at least square sense in the L2 sense, then the bad news is that you can't even get instance-optimality of order one unless you pick a number of measurements that's roughly the same as the length of the vector. Well, we don't want to do that. We're trying to have little and small. We don't want little and big. So this says that instance-optimality in L2 doesn't seem to be a viable concept. Now what happens in LP or L3, for Q between 1 and 2 is you sort of mediate from one to the other. Near 1 you have almost the same as L1, and then as you move to Q equal 2, measuring the error in LQ norms, the situation deteriorates until when you get to L2, it completely falls apart. Well there's some recovery. Although I said that instance-optimality in L2 is not possible, there is a little bit of recovery. And I think this is an interesting fact that you can recover some. I want to change the game a little bit. Well the game I've played so far has been you create once and for all a matrix phi, you create once and for all a decoder, and I use that system over and over and over again. Now I want to change the game. What I want to allow is a different game. That I give it, I want to think that I have a bunch of encoding schemes. For example, I have stochastic matrices, and I can grab any one of them that I want. And I start with a vector x, and I roll the dice and grab one of these stochastic matrices, and I use it for my encoding. And then I decode it in some way. And now I ask, okay, I know that I'm not going to be able to do for every x with certainty get instance-optimality, but maybe I can do it with high probability. And this turns out to be the case. And to describe this, I need to tell you what types of matrices, or what do you need about the stochastic matrices to make this go. So suppose I have a collection of random matrices. I want this collection to have two properties. The first is I want it to have the restricted isometry property of order k with high probability. Well we know that. We know how to create such matrices. That when we make a draw, the resulting matrix will have a restricted isometry property with high probability. So I want this. And the second thing I want is I want it to be, remember the restricted isometry property applied to vectors in sigma k. Here I want to have a bound in this property with high probability. Namely that for every vector x in Rn, I want with high probability that a random draw of this matrix when I look at phi omega applied to x, its norm is controlled by, let's say, two times the norm of x. Now we know in general, given an x, in general we can't, or let's say given the matrix phi of omega, when I start putting in different x's, for some x's this is going to be big. The norm is going to be big. But we also know that the probability that it's going to be big is small. Remember those concentration and measure inequalities? They said the probability of this was small. So it turns out you can have matrices that satisfy this property, the Gaussian, the Bernoulli. They all satisfy these two properties. And now here's the theorem. The theorem says suppose you have a collection of matrices and they satisfy those two properties, then with high probability you're going to have instance optimality. Namely take your x, draw your matrix at random, apply your encoding and decoding, then with high probability you're going to have this estimate. So should we be happy with results in high probability? I think so. I mean, if you know that you can do something with probability 10 to the minus 10, that's probably more certain than anything else you'll ever do in your life. So I think results in high probability are interesting. Okay. And the range of k, again, is the same range as the best that we could expect. An unfortunate part of this is that this presently stands, the way we prove this theorem is that we use this least squares minimization, this thing that we can't really do. I mean, it's this combinatorial search over all these spaces. But I believe that we probably can use the L1 minimization. We just haven't proven that. Okay. I'm done with the content of the talk, but I want to make now some comments about, okay, where are we? What are we trying to do? Where are we? So I'm going to talk about three subjects here, amplify a little bit, sort of to give you an idea of where we sit. So let's talk about optimum matrices. My whole thing was to what's optimality and to build optimal matrices. The optimal matrices that we've built, whenever we ask for optimality in terms of performance of sparse vectors and stability in addition, the only way we've known to do this is through probabilistic methods. So it's interesting to ask, could we do this through deterministic constructions? Now I want to be a little careful on this because this needs some thought of what do you mean by deterministic? Because let's think about our Bernoulli matrices. We know we can create matrices with ones and zeros that are optimal. We prove that by a probability technique. What we need is to have this restricted isometry property. On a matrix, I can check the restricted isometry property, not realistically check it. It would take forever to check it, but if you say deterministic, yeah, we had enough people working all the time, we could check it. So I could start writing down all zero one matrices, there are only a finite number on them, check for restricted isometry property until I get one and I've got a matrix that works so I have a deterministic construction. But this isn't fair. To me, this isn't a deterministic construction. So I have to be careful about what I mean by deterministic construction. And maybe what we need to do is make this a precise statement, what we mean by deterministic. Maybe it means that you have an algorithm that in polynomial time will assemble a matrix that has the properties you want. We at present, we have no deterministic construction that gives the type of performance that we want, that gives the high performance. That is, we have to take very few samples. Remember we had to, we could get by our, if the sparsity was K, then the number of samples N was just a little bit bigger than K. What you can do, and this is very poor man proof, I mean it's very easy to do this, is you could prove that through using matrices from coding or finite fields, you could prove that you could construct matrices, really construct them, that would give sparsity performance for K less than square root of N. Well that's far from N. So we're far, a long ways away from having good performance there. Okay, second topic I want to say a few words about is computational issues. I have tried, I have made the following, there were three ingredients that came up when we started discussing optimality. Second was we wanted optimality in terms of the sparsity level, we wanted to be able to handle as sparsity signals as possible. Okay, a second one that came up was the number of computations you needed once you implemented your algorithm. And the third one was stability. And these things play against one another. When you try to push the sparsity level up, then you have problems with computations and you have problems with stability. And it has not been sorted out yet exactly what the story is here. What are the best results? So if you, if you fix the sparsity level, are you forced to have spend so much time decoding? The best we know right now is L1 minimization decoding. This is linear programming, could be order N cubed. Can we do better? This is not understood. There are things that are being done and unfortunately I didn't have time to mention it in theoretical computer science. Gilbert Mutuchrist, Nann Strousen, many other people. I mean, this is a problem that really fits well in theoretical computer science. I'm sure if there's any people from that area in the audience who are saying, well, we do this all the time. And, and it's sort of true. And they have results, but they don't come close to answering this triad of questions, how these play off. They introduce other aspects of computation into the picture. Finally, we're really interested in building practical schemes, right? I mean, we're not just doing this for, for the fun of it, although it's nice interesting mathematics. We would like to have analog signals. Our original problem was we have these cell phone conversations or our, our radar chirps out there and we want to be able to find them and grab them. That's the whole purpose of compressed sensing. So we really want analog, we want to handle analog signals. So how do we go from discrete to analog? This is not completely clear. The, the waters are very muddy here and it's not sure that we're going to really be able to do this in the way that we would like. And moreover, I want to say that even if you mathematically could conceive of a system that would work in analog, you have to keep in mind that we want to build this system. So we're going to put it into hardware. We're going to, when, when you start putting things into a circuit, you know, all of a sudden it's not just like even numerical computation. You have more serious problems to deal with. So to get the big bang that we're hoping for in signal processing is not there yet. But I do want to mention that there are already impressive applications of compressed sensing, in particular, Emmanuel Candez in tomography. And I'd like to recommend that you go to Emmanuel's talk that is taking place here, which takes place last Friday. Okay. One last thought, if you're at all interested in compressed sensing or these ideas, there is a website at Rice in the Electrical and Computing Engineering and it has fairly up to date all the papers dealing with compressed sensing. There's surely more than 100 of them now and it's an interesting area and if you want to, you can easily get into it from that perspective. Thank you very much for your attention.
 It gives me great pleasure and it's an honor to introduce Ron DeVore, who over the last  few decades has been pioneering and introducing a variety of methods of functional approximation  and approximation theory in all of its ramification.  Ron's contribution covers a broad range of issues and problems, deep analytical issues  related to methodology for functional approximation.  You will hear sort of the outcome and import of this kind of work in his talk later, but  let me just say that in our day functional approximation is sort of the key to the digital  age in some sense.  Machine learning and functional regression, it's sort of the main straightforward elementary  question that everybody asks, but nobody knows how to answer, and hopefully we'll learn something  in the next few minutes.  Later aside, I just learned a few minutes ago that we are actually mathematically related.  He's my nephew, mathematical nephew.  We graduated under Bojanic, who is my mathematical brother.  We are both students.  We were both students of Karamata.  So the world is small in mathematics and beautiful.  So thank you, Uncle Ralfi.  Ralfi.  I think we're all aware that most scientific problems can't be solved exactly, so we have  to resort to some sort of approximation or numerical computation.  So we create algorithms to try to resolve these problems, and my talk is concerned with is  there a way to tell whether an algorithm is optimal or near optimal?  This would be not only a wonderful result mathematically, but would lead us in the right  directions of how to proceed computationally.  So to answer that question, of course we're going to have to describe exactly what we  mean by optimal.  But I want to put another caveat into this.  You may think of an algorithm that sounds very good to you in theory, but you have to  remember that whenever we implement such an algorithm, we do it on a computer, and we're  subject to machine errors and even errors in entering the data.  So it's not just a question of designing an algorithm, but the algorithm should be numerically  stable, and we'll see that ingredient coming up in my talk.  So there are many settings that I could describe this, and in fact with my colleagues, especially  with Albert Cohen and Wolfgang Dahman, we've studied this problem of optimal algorithms  in many different settings.  But after some thought about what I should do here in this talk, I've decided to speak  about the, in the setting, this clicker is not working very well, so I think I'll stand  here maybe to do better.  I decided to concentrate on an area which is called compressed sensing.  This is a relatively new area in mathematics.  It's getting a lot of interest these days.  And the motivation for this is the problem of sensing real world signals.  So the usual paradigm for signals is that they're band limited functions, and the way  to sample them is through what's called Shannon sampling theorem or Nyquist sampling.  The unfortunate state of affairs is that many of the real world signals that we're interested  in are broad banded. This means that their Fourier transform has large support, and that  this type of sampling cannot be implemented in practice.  So compressed sensing is looking for a way around this to find another way, an innovative  way of sampling that defeats the roadblock of Shannon sampling, and samples signals  more at the information content in the signal rather than at the, determined by its bandwidth.  So my reason for choosing this subject is there are three reasons.  One is that it's very audience friendly as you'll see.  It'll be easy for you to grasp the main ideas that are here because it's basically going  to be an exercise in finite dimensional geometry and linear algebra.  The second reason is that the subject itself interfaces a lot of areas of mathematics,  most notably probability theory, theoretical computer science, and functional analysis.  But the real reason is that by discussing this subject, I'll be able to introduce the  types of optimality that are important in numerical computation and all the nuances.  As we'll see, I will introduce three types of optimality when we discuss this problem,  and it's sort of the full range that I know of that occurs in numerical problems.  Okay, so let's begin.  Well, really in compressed sensing, we're interested in analog signals, but the subject  is fully developed only for discrete signals.  So I'm going to stick to discrete signals until the very end when I make a few comments.  So our problem is the following.  We have a vector x and rn, and here capital N is large.  So for example, if we were dealing with images, n would be the number of pixels, it would  be of the size one million, or even more.  And the problem is that the game we're going to play is we're allowed to ask small n questions  about this signal.  And these should be non-adaptive, by which I mean I don't want you to have that the  one hundredth question depended on the answers to the previous 99.  You must set the questions down in advance.  I'm even more restrictive because what I mean by a question is going to be very specific.  A question to me is that it's going to be the inner product of this vector x with some  other vector.  So you can choose the vector v, and that will give you then by taking the inner product  of v with x, you'll get the answer to your question.  So question is very specific in this case.  It's linear information.  Now, we can always represent this kind of a system by an n by n matrix phi.  The rows of this matrix would be these vectors v that we're taking the inner product with  respect to.  And capital N is the length of the vector x.  So typically, capital N is very large, and small n is very small.  So we have short and fat matrices phi.  So you see in this world, the Brazilian diets, this is not applied.  We're interested in short, fat things.  Okay.  My question is, can we say what are the best matrices to use in such a system?  Or if not best, near best, or good matrices?  Now here, good is going to mean what?  It means that once we extracted this information y, which is phi of x, so x appeared before  us.  We asked these questions.  We get this vector y, and now x disappears from our site, and all we can use is y.  And what we're interested in is how well can we recover x from this information y?  How well can we approximate x?  This is roughly speaking what we mean by good.  Now there's going to be two aspects to this problem.  One is, does y contain enough information about x?  And the second is, how do we extract this information from y?  That part is called the decoding, and both are serious issues.  Okay so let's think about what the problem really is here.  So we have matrix phi that's mapping a large dimensional space, Rn, into a small dimensional  space.  So naturally, a lot of vectors are mapped into the same image.  In particular, if you look at all the vectors mapped into zero, the null space, this is  a large space.  It will have dimension capital N minus little n.  Since capital N is big and little n is small, this is a large dimensional space.  We have a lot of collapsing going on.  If you took any vector y and asked which vectors are mapped into y, they can be described  as just a hyperplane where you take any one of the vectors mapped into y and then add  the null space to it.  This gives you a hyperplane, and this entire hyperplane has the same information y.  So these hyperplanes which are denoted here in yellow, so this is a given hyperplane for  a given y, all the points on this hyperplane are mapped into the same information vector  y.  Right, so there's tremendous collapsing here going on.  Now a decoder, what is it going to do?  It's going to take y and it's going to map us back into the big space, rn.  So after it does that, we get our approximation x bar to our original x.  Now at this point, you should be very pessimistic that we're going to do something here because  of this collapsing.  So this vector x bar would have to serve as an approximation for all the vectors in this  big hyperplane that was mapped into y because none of those are distinguished in any way  by this encoding decoding.  They're all treated the same way.  Well indeed if we had to work in complete generality, we wouldn't be able to say much  about this.  However, our understanding of the world of real signals is that they have some structure  and we want to utilize this structure or take advantage of this structure in going forward.  So our paradigm or the way we think of real world signals is that they're sparse in some  sense, that there's a set of building blocks and if you represent this signal with respect  to this set of building blocks, it either has a very small representation or it can  be approximated well by a few terms.  Now building blocks, in my case, I'm going to restrict to the case where the building  blocks are a basis for our end.  You could think of this in other settings as well.  In some problems, we may know the basis in which the vector is sparse.  In other problems, we might not know the basis.  So this is an important consideration.  I'm going to actually proceed for the moment assuming that the basis is known to us and  then later I'll say something about what we do when the basis is not known.  So here's going to be my first way of measuring optimality.  So I'm interested in telling you that when we tackle a problem like this, let's try to  make a precise mathematical problem to see what we mean by a best or near best algorithm.  So here's my first notion of optimality.  I want to deal, I'm going to assume that my signal x is sparse.  And I'm going to, for now, I can assume that the basis in which it sparses a canonical  basis.  This is not, I'm not pushing anything under the rug because you can always make a basis  transform to get to this case.  So what do I mean by a sparse vector?  I mean that it has few non-zero components.  By let's, given a vector x, the support is the number of i's, the indices for which the  vector is non-zero.  And I form this class, sigma k, which are all vectors whose support is less than or  equal to k.  So this is a non-linear set.  I mean it's not a linear space.  It's a union of hyperplanes.  You can take, for example, take any index at t of size k and look at the set of vectors  supported on t.  And then you union these.  So you can visualize this as you have a lot of hyperplanes and the union of these is the  set sigma k.  So I'm going to assume my signal is in sigma k.  And I'm going to ask, given that information, how can you best design an encoding-decoding  scheme?  So here's the problem.  No one is fixed, little k is fixed.  And I ask, what's the smallest n such that I can build a matrix phi and a decoder that  will capture every one of those guys exactly?  And I'm going to tell you right off hand what the answer is.  The answer is n equals 2k.  Now this is sort of morally justified if you think about it.  Why? Because a vector in sigma k is determined by 2k pieces of information.  You need to know the k positions where it's non-zero.  And then you need to know the value of x at these k positions.  That's 2k pieces of information.  So one would think you can do this with 2k samples.  But of course, my definition of question was very specific.  I have to take inner products.  So it's not absolutely clear that I can do this.  But as we will see, we can.  So to describe what matrices, so my goal now is to describe to you what are the good matrices  for this problem?  What matrices will work?  So I want to think of my matrix as built from its columns.  So this matrix has columns v1 through v capital N. And I'm going to introduce a property  of the matrix that if the matrix has this property, it will allow a solution to this  problem.  And to describe this property, what I look at is sub matrices of phi.  And the sub matrices I get is I take any t columns, any set t of indices, column indices,  and I look at the sub matrix gotten from phi by just keeping those columns and throwing  away everything else.  I'll denote this by phi t.  And I will use this matrix frequently.  And associated to that matrix is this Gramian matrix, phi t transpose phi t, which is the  matrix of inner products.  This matrix is a square matrix.  It's symmetric and it's non-negative definite.  And now here's the theorem.  So suppose you have an N. And you want to check whether a matrix of this size will  capture every k sparse vector exactly.  I claim this will be the case if and only if one of these four conditions hold and they're  all equivalent.  So the first condition is that the intersection of sigma k of these three conditions, since  this is the requirement, that we reproduce all signals exactly.  So the first condition is that sigma 2k intersects in all spaces zero.  Now notice I'm trying to recover k sparse signals, but here the condition is 2k.  Now why do we have that?  It's very intuitive.  We need to avoid the case where two signals of length k are both mapped to the same y,  because then we couldn't find which one it was.  While two signals are mapped to the same y, I can take their difference and that'll be  mapped to zero.  Well that's something in an all space, but their difference will typically be a signal  of length 2k.  So that's how I get to this condition.  And equivalent conditions are the following and these are the simple linear algebra.  Another way to state this is to say that these matrices phi t that I form, no matter what  selection t of columns I pick, as long as I pick 2k columns, that this matrix should  have rank 2k.  Another way to say it is that this Gramian matrix should be invertible.  So if I can find a matrix phi with these properties, then I know it will recover two k sparse signals  exactly.  So the question becomes, can I construct such matrices?  Can I write some down?  Maybe I'm stating something you can't find such matrices.  Well what do we need to do?  We need to create vectors in r2k.  I'm claiming I can do this for n equals 2k, right?  But I need to be able to create a lot of vectors because capital N is arbitrary.  It can be arbitrarily large.  I need to create a lot of vectors in r2k and I need that whenever I take 2k of them, they're  more or less independent.  They have to be linearly independent.  Well here's a quick and easy solution to this.  Take the Van der Mann determinant.  Take distinct points x1 through xn and form the matrix where in each column you have the  successive powers of the given xj's.  Then we know that if we pick a sub-matrix consisting of 2k columns of this, we get a  square matrix which corresponds to polynomial interpolation.  This matrix is invertible.  So this matrix satisfies the properties I wanted.  So it solves the problem.  But I haven't discussed how I'm going to do decoding.  Now you can think of many ways to do decoding.  I'm going to just put down a naive decoder.  A naive decoder would be I'm looking, now I'd like my decoder to be defined for every  y, not just for the y's that are images of k-sparks vectors because in practice I won't  have that information.  So I want my decoder to be defined for every y.  So here's something you could do.  Given y, you could look over this set sigma k and try to find a z that best fits y, that  is that y minus phi of z is small as possible.  This is a least squares problem of trying to find a least squares fit to y from the image  of sigma k under your map phi.  You can decouple this problem.  You may look at these spaces that we have introduced before where t is a set of column  indices of size k and we solve the least squares problem on this set t.  And of course we know we can actually write down what the solution is.  This is more Penrose inverse where we have to apply this Gramian matrix, the inverse  of it, to phi ty.  So this already shows you where we need somehow that this matrix is invertible.  We need non-singularity.  And then we just search over all these t and find the t star that gives us the smallest  residual.  In the case that our vector was the image of a k-sparks vector, there will be a t for  which this is zero.  And it will be a unique t and we'll pick up that t by this decoding process.  Okay, sounds great.  We solved our first problem.  Can we now go home and have a glass of cognac and feel happy that we've done something positive?  Well, I'm sorry to tell you that there's bad news here.  None of us will be alive when this decoding is finished.  Okay?  So this decoding that I told you, this naive decoding is horrible.  If you had a problem where n is like a million, like an image processing and little n is like  a thousand, you can imagine you're going to have to take capital n things taken little  n at a time to check.  But a more serious problem, it turns out, is that the decoding is not stable.  Now you can actually fix the first problem.  You could find another matrix.  For example, you could take the discrete Fourier matrix, first n rows of a discrete  Fourier matrix, and then by some clever decoding, you can do this whole problem in this number  of operations, which is a reasonable number of operations for us.  So this isn't the real problem that the decoding takes forever.  But the fact that the decoding is unstable is a real problem.  And let me just explain this to you that we're really stuck here in some sense, because  suppose I'm a generous man.  Suppose I tell you where the vector x, what its support is.  Suppose I tell you, OK, it's supported on t.  Now you only have to find, you don't have to find these indices of support.  You only have to find the vector x.  How would you do this?  Well, you have to apply this inverse operator.  And this, of course, we know in the examples that I created that this is very poorly behaved.  And in fact, now, if I stick to wanting to have only 2k samples, I'm stuck.  I cannot do this in a stable way.  OK, so what seemed like an easy mathematical solution, a lesson to us, is what seems mathematically  OK may not be numerically reasonable.  Well this brings us to actually where compress sensing began.  I'm not telling you this in chronological order.  Compress sensing began about three or four years ago in the work of Immanuel Kandez,  Terry Tao, sometimes Justin Romberg, and in a parallel development by Dave Donohal.  And what they do is they show how to create matrices that not only will recover 2k sparse  vectors, but do it stably.  To do this, one has to pay a little bit.  One has to increase the number of questions, and we'll see how that plays out.  So they make, and I'm going to talk about the Kandez-Tao development, which I think  is very elegant and very much to the point, they make two important contributions.  First, they're going to tell us what are the good matrices, or they're going to tell us  a sufficient condition that a matrix be good.  And secondly, they're going to handle this problem of decoding.  They're going to show us how we can do fast decoding.  So how does this play out?  So the property that they're going to introduce for the matrices is the following.  Remember these matrices, phi t and these Gramian matrices, they're going to require, we say  the matrix satisfies a restricted isometry property of order k.  If whenever we pick k column indices, the eigenvalues of this matrix are in 1 minus  delta 1 plus delta, and here delta is between 0 and 1.  So the important thing is the eigenvalues are staying away from 0, and they're staying  bounded.  So this Gramian matrix now is not only invertible, but it's nicely conditioned.  So that's the new ingredient.  You could restate this property as follows.  You could just say that what you need is that whenever you apply phi to any vector in sigma  k, the resulting vector has norm roughly the same as the norm of x.  That's where the word restricted isometry comes in.  If we could take delta in this equal to 0 in this formulation, then it would mean that  we have an isometry.  But we give up a little bit, and we have to give up a little bit.  Now what about the decoding?  They suggest for decoding that once use L1 minimization in this way.  So suppose we have our vector y.  To decode what I'm going to do is look at all x which map into y, and then I'm going  to pick the one which has minimum L1 norm.  This has a very simple geometrical description.  If I think of y, it determines this hyperplane, and this is the ball in little L1.  I take this ball very small about the origin, and I blow it up like a balloon until it hits  this hyperplane.  When it hits this hyperplane, that would be the solution.  And generically, it usually happens that it hits at a unique point.  Now here's the main result.  By the way, an important part of L1 minimization is it can be solved by linear programming.  So we have a modest investment in the decoding.  For example, we always know we can decode with order capital and cubed operations.  And here's the main result of Candidaz and Tal.  They say if you have a matrix phi that satisfies this restricted isometry property, now for  3k, before I had 2k, probably you could replace 3k by 2k and everything if one studied this  enough.  But once you have this property, then you can find a matrix with this.  If you have a matrix with this property, then you use L1 decoding, you'll capture every  x exactly.  That's in sigma k.  Every sparse guy will be captured exactly, but moreover, the decoding is stable.  So now our problem is, OK, we have a new matrix problem.  Then we find matrices phi that satisfy this restricted isometry property.  And how do we do that?  So how do we build good matrices?  Well, I claim that for given n and n, we construct such matrices.  So we'll satisfy restricted isometry property of order k.  For every k less than n divided by log n.  Now all of a sudden you see a log n appear in this, which hadn't appeared before.  This is the price we pay for stability.  You can actually show that you can't discard this log n.  It must be there.  If you want this restricted isometry property, you're stuck with this log n.  So we give up a little.  Before we could take k equal to n over 2.  Now we can't quite take k at large, but only sacrifice by this logarithm.  Now, the question I want to discuss next in the next few minutes is, how can we construct  such a phi?  What do we want to do?  Well, we're back to our problem of creating a matrix with a lot of columns, matrix and  our vectors and our n, and a lot of them.  Now we want them not only to just be that whenever we pick 2k of them or k of them that  they're linearly independent, we want them to be very far from being dependent.  We want them to be almost orthogonal in some sense.  I'm going to mention to you now three ways you can do this.  All of these constructions are probabilistic.  The first way is we could look at the unit sphere and our n and to choose our columns  we simply pick at random according to the uniform distribution.  Independently at random we pick our columns.  This would work.  Another way is we could create a stochastic matrix.  What we could do is take one of your favorite distributions, let's say the Gaussian distribution  with mean zero and variance one over n.  And for the first entry just make a draw of this distribution and whatever number you  get stick it in as your first entry in the matrix.  And then independently repeat this every time you want an entry and thereby fill out your  matrix.  You could do the same thing with a Bernoulli distribution with plus or minus ones or zero  and ones.  Only one of these constructions will work and will give you a matrix phi with high probability  it will work.  That is you do this realization then you can prove that with overwhelming probability  the resulting matrix will satisfy the RIP property for this range that I advertise and  I tell you this is the best range you possibly can have.  So we have a way to construct these matrices.  I want to make one point very clear here where probability is coming in because sometimes  people get confused on this point.  We use probability to construct this matrix phi to prove the existence of phi.  What we know is we do this construction and almost every matrix we get is going to have  the property we want.  Once we have a matrix in our hands the algorithm is constructed.  There's nothing probabilistic in the algorithm.  We apply phi, we apply L1 minimization and so on.  Probability is only used to construct this.  Because this RIP property is so important, it's of interest to have nice methods to verify  RIP.  Typically one would think that this is an eigenvalue problem and I'll approach it through  eigenvalues.  What I want to point out here is that actually there's an elementary track to verify the  RIP property in the settings that I just described, for example in the Gaussian or Bernoulli  case.  So suppose I generate a stochastic matrix by taking draws of some random variable R and  fill out the matrix in this way.  Then the first thing to observe is you can just write out what phi omega of X is and  expand it and take the norm and you'll see that the expectation of this vector is norm  of X squared.  So I'm going to try to draw out what distributions, what probability distributions would work  in creating these matrices that satisfy the RIP property.  So this is always going to hold for any distribution, but we're going to need more.  What we need to have more is that when we look at a given X and we look at phi omega  of X that it's with high probability close to this mean.  So we need a concentration of measure inequality, namely that for any given X we have this estimate  on the probability that these two things differ by a delta times the norm of X. This concentration  of measure inequality will be satisfied for Gaussian, Bernoulli, many distributions and  it's usually an easy matter to show that this property is satisfied.  And I claim that whenever you have this, then the resulting matrix that you create from  this random process will satisfy the RIP property for the range of K that you want.  And maybe it's not so important that this is true, but what I want to emphasize is the  proof because there are other ways, as I mentioned through Eigenvalue analysis, to prove this  result, but there's actually this has a very simple analytic proof that I want to just  quickly draw upon.  So how would I prove this?  What I do is, see I need to verify on sigma K, I need to verify this restricted isometry  property which meant that the norm of phi applied to X is like the norm of X for all  X and sigma K. That's what I need.  So what I do is first I cover the unit sphere in sigma K to some tolerance, let's say delta  over 4.  Then I know that with this covering, this net of points, for each point in the net I know  this, this is a concentration of measure inequality, told me this.  With high probability this is true.  And since it's with high probability I can apply a union bound and get that as uniformly  true for all Q and P. And of course you've got to be careful of how many points do you  throw in there, but do this.  And then finally, you need to extend it to every vector in sigma K, not just the points  in this cover P. And to do that is quite trivial.  For example, the upper estimate, you can estimate norm of phi by saying, well I'll approximate  X by a Q from my net.  And then I have this trivial estimate for norm of phi of X. And the second one is under  control because of this.  And the first one, well any matrix on a finite dimensional space is bounded, so there is  some bound.  But of course I won't get the best result this way with M. But then I can bootstrap  this.  Because once I get a bound for M I can repeat the process.  And if you do this you'll arrive at this theorem.  Okay, now I talked in the beginning about, you know, the basis is known to me.  What can we do if the basis is not known?  And people in this subject talk frequently about this concept of universality.  Well, these random things, they work for any basis.  But I caution you to be a little bit careful on this point.  Let's make this precise.  To make this precise, if I start with any collection of bases, and if the number of  bases in this collection is controlled by e to the cn, n is this, you know, remember  the number of rows in a matrix.  Then this same argument I gave probability argument can be used to show that you can  prove the existence of one matrix phi that will satisfy RIP with respect to every one  of these bases simultaneously for the hope, the range of k that you would like, the biggest  range of k.  The problem, so this is a sense that you can say universality.  You can say that if you draw at random a matrix with high, and if you have these bases known  in advance to you, this collection of bases, and with high probability, this will uniformly  handle and encode well, spar signals with respect to any one of these bases.  The problem is that we would need to know the basis to decode.  This is good for encoding, but not good for decoding.  Well, some people talk about, well, I'll encode and then 20 years later I'll decode.  You know, I'll put it in a time capsule.  I don't know too many applications in which people want to wait 20 years, but.  Okay, that was my first example of optimality, and it was the simplest.  Now I want to go to more general.  So we know that real world signals aren't typically going to be spars in the sense that  I mentioned to you.  They're not going to have support that's only a finite number of terms or a small number  of terms.  They're typically going to have a lot of entries in them.  But in general, most of these entries are going to be small.  And our idea, remember as we began, was that we said that we view signals as being approximated  well by spar signals.  That is, they have only a few components that are essential to capture them, and we only  need to find those components.  So to make this meaningful discussion, I want to talk about how a signal X can be approximated  by spar signals.  So I introduce this measure here of error.  So you can choose whatever way you want to measure error.  The norm X maybe for this talk just take X is the least squares L2 and and then given  a vector X, we look at how well it can be approximated by vectors which have support  K. So a vector is nice if there's a vector of length 100, support 100 that approximates  it well.  And we think real world signals are nice.  So mathematically, we're interested in classes of signals in which we know that they can  be approximated well.  And typical way to get such classes are the following.  You could take the unit ball of an LPN.  So for example, if you want to measure error in LQN, and I suggest you to take Q equals  2 here, then any vector in LPN will be approximated by K terms to this accuracy.  Now look at this a minute.  As P gets small, now here small P is better than large P. Small P means that the vector  X, the components of it die off very fast if you order them in size.  And you'll have this decay, K to the 1 over Q minus 1 over P. So this can be very, very  fast if P is very small.  For example, if Q is 2 and you take a vector in a unit ball of L1, then it can be approximated  by K terms to accuracy 1 over square root of K. This is easy to see and I don't stop  to prove that for you.  So what would I mean by optimality in the second sense?  Optimality now will mean I have a class of signals, capital K, and I have some norm in  which I'm going to measure distortion or error.  And I say that the sensing system is optimal if it performs like this min-max.  So this min-max says, OK, let's look at how well we can possibly perform.  This is the best performance we can make by systems of size n.  That is, the matrix has a dimension, the little n by capital N.  We take such a matrix and a decoder and we apply this to a given X and we achieve a certain  error.  And now we look at the super, the errors over all X in the class K. So this is a usual min-max  way of doing things.  We look at the worst error and then we look at the system, the phi and delta, that give  us the best worst error.  So we take the infimum over all these.  And now I would say a system is near optimal if I have this kind of an inequality with  some absolute constant C0.  If I had C0 equal 1, I would have optimality.  There would be nothing better, but usually you have to give up a little something.  So this is interesting to understand for a given class of signals, how does this en behave?  What's the target?  What's the best we could possibly do?  That's the usefulness of en.  It tells us the best we could do.  And the story is that this, the behavior of enk for classical sets K in our n is known.  It was not only known, but it was done a long time ago in the 1970s and early 80s in the  context of finite dimensional geometry or geometry of Banach spaces.  The main players were Gaussian and Gluck skin.  And I just mentioned one result.  For example, if I look at the unit ball in L1n and ask what's the best encoding decoding  I can do with that in L2 when I measure the error in L2n, then here's the estimates.  Notice the estimate from above and below are the same order.  The only thing that differs are these constants, C0 and C1.  The estimate is basically 1 over square root of n.  1 over square root of n would be what we could do if we chose the n best terms of the vector.  And we pay a little price here, this log of N over little n again.  Remember we had that before.  And this price we know is necessary to pay it because we have a lower estimate.  Now I'd like to say so why was this not a solution to the compressed sensing problem  back in the 70s?  The thing that these players did not look at was practical decoding or how can you implement  this encoder?  In fact, they never touched the problem of decoding.  And a contribution of Candidaz and Tao and Donohull is that they show that if you have  a matrix with the restricted isometry property and you add L1 minimization, then it will  give near optimal performance.  It'll achieve this type of performance.  Now not for every Q and P, but for certain Q and P. For example, if you want to measure  error in L2, then for every P less than or equal to 1 it'll be optimal.  So they give a practical scheme and once you have this scheme, now you're happy because  you know that you're doing the best you possibly can.  Let me say just a couple words about this because there's, you know, I wouldn't say  a controversy but I talk to people in functional analysis and they say, look, these guys in  compressed sensing, they're just doing our work over again and so on.  This is not a fair statement because as I mentioned to you, they did not discuss the  practicalities of the decoding.  So their results are useful and interesting but didn't discuss practical decoding.  What is the most useful about the results is the lower bound, it turns out, because  it tells us we can't do any better.  But it doesn't give us real practical schemes to do the encoding and decoding.  So I think both camps have made a great contribution.  Now I want to go to my third and last way of measuring optimality.  And it, to me, if you were in the States, you would call this the Cadillac of all ways  of measuring optimality.  And Germany would be, I guess, the Mercedes Benz.  I don't know what it is in Spain.  But anyway, what I want to have optimality in a very strong sense, here's the sense.  I'll say that this encoding-decoding scheme is instance-optimal of order K. If the following  holds it, every vector I choose.  I don't need it in any class or anything like that.  Take any vector.  When I encode and decode, it performs like K term approximation.  Well, this would be terrific if you can do this for large values of K. For example, if  you have this and you can return to your first problem, and if you take a vector which is  K sparse, this will be zero.  So you'll encode, you'll get it exactly.  It will also include all the results on optimality for classes.  So this to me is like, you know, if you could get this, it's like a dream, but if you could  get this, this would be the best possible result.  And the problem is, let's suppose you fix the size of the problem, the size of the system  you're going to use, the number of rows you're going to use in the matrix.  Then the question is, how big can you take K and achieve this instance-optimality?  And with L Bear Cohen and Wolfgang Dahman, we've solved this problem for all, if you  want to measure the error in any LQ space as long as Q is between one and two.  And I'm just going to mention a couple of these results.  So there are two things to say.  One, the good news.  The good news is that you have instance-optimality in L1, and you have it for K of the range  that you're the most you could expect, namely K less and log n.  Moreover you can build this system and use L1 minimization as a decoder so you can actually  build a stable system and have this instance-optimality.  And this result, although we prove it in our paper, you can, it essentially follows from  the work of Candez and Tao.  Now that was the good news.  Now here's the bad news.  The bad news is that if you ask for instance-optimality in L2, you want to measure error in at least  square sense in the L2 sense, then the bad news is that you can't even get instance-optimality  of order one unless you pick a number of measurements that's roughly the same as the length of the  vector.  Well, we don't want to do that.  We're trying to have little and small.  We don't want little and big.  So this says that instance-optimality in L2 doesn't seem to be a viable concept.  Now what happens in LP or L3, for Q between 1 and 2 is you sort of mediate from one to  the other.  Near 1 you have almost the same as L1, and then as you move to Q equal 2, measuring the  error in LQ norms, the situation deteriorates until when you get to L2, it completely falls  apart.  Well there's some recovery.  Although I said that instance-optimality in L2 is not possible, there is a little bit  of recovery.  And I think this is an interesting fact that you can recover some.  I want to change the game a little bit.  Well the game I've played so far has been you create once and for all a matrix phi,  you create once and for all a decoder, and I use that system over and over and over again.  Now I want to change the game.  What I want to allow is a different game.  That I give it, I want to think that I have a bunch of encoding schemes.  For example, I have stochastic matrices, and I can grab any one of them that I want.  And I start with a vector x, and I roll the dice and grab one of these stochastic matrices,  and I use it for my encoding.  And then I decode it in some way.  And now I ask, okay, I know that I'm not going to be able to do for every x with certainty  get instance-optimality, but maybe I can do it with high probability.  And this turns out to be the case.  And to describe this, I need to tell you what types of matrices, or what do you need about  the stochastic matrices to make this go.  So suppose I have a collection of random matrices.  I want this collection to have two properties.  The first is I want it to have the restricted isometry property of order k with high probability.  Well we know that.  We know how to create such matrices.  That when we make a draw, the resulting matrix will have a restricted isometry property  with high probability.  So I want this.  And the second thing I want is I want it to be, remember the restricted isometry property  applied to vectors in sigma k.  Here I want to have a bound in this property with high probability.  Namely that for every vector x in Rn, I want with high probability that a random draw of  this matrix when I look at phi omega applied to x, its norm is controlled by, let's say,  two times the norm of x.  Now we know in general, given an x, in general we can't, or let's say given the matrix phi  of omega, when I start putting in different x's, for some x's this is going to be big.  The norm is going to be big.  But we also know that the probability that it's going to be big is small.  Remember those concentration and measure inequalities?  They said the probability of this was small.  So it turns out you can have matrices that satisfy this property, the Gaussian, the Bernoulli.  They all satisfy these two properties.  And now here's the theorem.  The theorem says suppose you have a collection of matrices and they satisfy those two properties,  then with high probability you're going to have instance optimality.  Namely take your x, draw your matrix at random, apply your encoding and decoding, then with  high probability you're going to have this estimate.  So should we be happy with results in high probability?  I think so.  I mean, if you know that you can do something with probability 10 to the minus 10, that's  probably more certain than anything else you'll ever do in your life.  So I think results in high probability are interesting.  Okay.  And the range of k, again, is the same range as the best that we could expect.  An unfortunate part of this is that this presently stands, the way we prove this theorem is that  we use this least squares minimization, this thing that we can't really do.  I mean, it's this combinatorial search over all these spaces.  But I believe that we probably can use the L1 minimization.  We just haven't proven that.  Okay.  I'm done with the content of the talk, but I want to make now some comments about, okay,  where are we?  What are we trying to do?  Where are we?  So I'm going to talk about three subjects here, amplify a little bit, sort of to give  you an idea of where we sit.  So let's talk about optimum matrices.  My whole thing was to what's optimality and to build optimal matrices.  The optimal matrices that we've built, whenever we ask for optimality in terms of performance  of sparse vectors and stability in addition, the only way we've known to do this is through  probabilistic methods.  So it's interesting to ask, could we do this through deterministic constructions?  Now I want to be a little careful on this because this needs some thought of what do  you mean by deterministic?  Because let's think about our Bernoulli matrices.  We know we can create matrices with ones and zeros that are optimal.  We prove that by a probability technique.  What we need is to have this restricted isometry property.  On a matrix, I can check the restricted isometry property, not realistically check it.  It would take forever to check it, but if you say deterministic, yeah, we had enough  people working all the time, we could check it.  So I could start writing down all zero one matrices, there are only a finite number on  them, check for restricted isometry property until I get one and I've got a matrix that  works so I have a deterministic construction.  But this isn't fair.  To me, this isn't a deterministic construction.  So I have to be careful about what I mean by deterministic construction.  And maybe what we need to do is make this a precise statement, what we mean by deterministic.  Maybe it means that you have an algorithm that in polynomial time will assemble a matrix  that has the properties you want.  We at present, we have no deterministic construction that gives the type of performance that we  want, that gives the high performance.  That is, we have to take very few samples.  Remember we had to, we could get by our, if the sparsity was K, then the number of samples  N was just a little bit bigger than K. What you can do, and this is very poor man proof,  I mean it's very easy to do this, is you could prove that through using matrices from coding  or finite fields, you could prove that you could construct matrices, really construct  them, that would give sparsity performance for K less than square root of N. Well that's  far from N. So we're far, a long ways away from having good performance there.  Okay, second topic I want to say a few words about is computational issues.  I have tried, I have made the following, there were three ingredients that came up when we  started discussing optimality.  Second was we wanted optimality in terms of the sparsity level, we wanted to be able to  handle as sparsity signals as possible.  Okay, a second one that came up was the number of computations you needed once you implemented  your algorithm.  And the third one was stability.  And these things play against one another.  When you try to push the sparsity level up, then you have problems with computations and  you have problems with stability.  And it has not been sorted out yet exactly what the story is here.  What are the best results?  So if you, if you fix the sparsity level, are you forced to have spend so much time  decoding?  The best we know right now is L1 minimization decoding.  This is linear programming, could be order N cubed.  Can we do better?  This is not understood.  There are things that are being done and unfortunately I didn't have time to mention it in theoretical  computer science.  Gilbert Mutuchrist, Nann Strousen, many other people.  I mean, this is a problem that really fits well in theoretical computer science.  I'm sure if there's any people from that area in the audience who are saying, well,  we do this all the time.  And, and it's sort of true.  And they have results, but they don't come close to answering this triad of questions,  how these play off.  They introduce other aspects of computation into the picture.  Finally, we're really interested in building practical schemes, right?  I mean, we're not just doing this for, for the fun of it, although it's nice interesting  mathematics.  We would like to have analog signals.  Our original problem was we have these cell phone conversations or our, our radar chirps  out there and we want to be able to find them and grab them.  That's the whole purpose of compressed sensing.  So we really want analog, we want to handle analog signals.  So how do we go from discrete to analog?  This is not completely clear.  The, the waters are very muddy here and it's not sure that we're going to really be able  to do this in the way that we would like.  And moreover, I want to say that even if you mathematically could conceive of a system that  would work in analog, you have to keep in mind that we want to build this system.  So we're going to put it into hardware.  We're going to, when, when you start putting things into a circuit, you know, all of a  sudden it's not just like even numerical computation.  You have more serious problems to deal with.  So to get the big bang that we're hoping for in signal processing is not there yet.  But I do want to mention that there are already impressive applications of compressed sensing,  in particular, Emmanuel Candez in tomography.  And I'd like to recommend that you go to Emmanuel's talk that is taking place here, which takes  place last Friday.  Okay.  One last thought, if you're at all interested in compressed sensing or these ideas, there  is a website at Rice in the Electrical and Computing Engineering and it has fairly up  to date all the papers dealing with compressed sensing.  There's surely more than 100 of them now and it's an interesting area and if you want to,  you can easily get into it from that perspective.  Thank you very much for your attention.
devore2006_0_6275.png devore2006_6276_9131.png devore2006_9132_10691.png devore2006_10692_10739.png devore2006_10740_12371.png devore2006_12372_12851.png devore2006_12852_16043.png devore2006_16044_16283.png devore2006_16284_16895.png devore2006_16896_16979.png devore2006_16980_26963.png devore2006_26964_27971.png devore2006_27972_29423.png devore2006_29424_29915.png devore2006_29916_32291.png devore2006_32292_33023.png devore2006_33024_33215.png devore2006_33216_33419.png devore2006_33420_33839.png devore2006_33840_34883.png devore2006_34884_38195.png devore2006_38196_40319.png devore2006_40320_40991.png devore2006_40992_43511.png devore2006_43512_45119.png devore2006_45120_45515.png devore2006_45516_48503.png devore2006_48504_49631.png devore2006_49632_54059.png devore2006_54060_56759.png devore2006_56760_58043.png devore2006_58044_59339.png devore2006_59340_61847.png devore2006_61848_66635.png devore2006_66636_66767.png devore2006_66768_66839.png devore2006_66840_68423.png devore2006_68424_72107.png devore2006_72108_74471.png devore2006_74472_75167.png devore2006_75168_76415.png devore2006_76416_76631.png devore2006_76632_77675.png devore2006_77676_78083.png devore2006_78084_78959.png devore2006_78960_81167.png devore2006_81168_81959.png devore2006_81960_82559.png devore2006_82560_84563.png devore2006_84564_84815.png devore2006_84816_85343.png devore2006_85344_85523.png devore2006_85524_85655.png devore2006_85656_86327.png devore2006_86328_87703.png