October 6.079/6.975 2930, 2009.
S. Boyd & P. Parrilo
Midterm exam
This is a 24 hour take-home midterm exam. Please turn it in to Professor Pablo Parrilo, (Stata Center), on Friday October 30, at 5PM (or before). You may use any books, notes, or computer programs (e.g., Matlab, CVX), but you may not discuss the exam with anyone until October 31, after everyone has taken the exam. The only exception is that you can ask us for clarification, via email. Please address your emails to both professors and the TA. Please make a copy of your exam before handing it in. When a problem involves computation you must give all of the following: a clear discussion and justification of exactly what you did, the Matlab source code that produces the result, and the final numerical results or plots. Matlab files containing problem data are available on Stellar. All problems have equal weight. Some are easier than they might appear at first glance. Be sure to check your email and the course web site on Stellar often during the exam, just in case we need to send out an important announcement. And one technical comment. For problems that require you to work out a numerical solution, you are welcome to use a solution method that involves solving more than just a single convex optimization problem. (Of course, only when this is necessary.)
1
1. 2D filter design. A symmetric convolution kernel with support {-(N - 1), . . . , N - 1}2 is characterized by N 2 coefficients hkl , k, l = 1, . . . , N.
These coefficients will be our variables. The corresponding 2D frequency response
(Fourier transform) H : R2 R is given by
H(1 , 2) =
k,l=1,...,N
hkl cos((k - 1)1 ) cos((l - 1)2 ),
where 1 and 2 are the frequency variables. Evidently we only need to specify H over
the region [0, ]2 , although it is often plotted over the region [-, ]2 . (It won't matter
in this problem, but we should mention that the coefficients hkl above are not exactly
the same as the impulse response coefficients of the filter.)
We will design a 2D filter (i.e., find the coefficients hkl ) to satisfy H(0, 0) = 1 and to
minimize the maximum response R in the rejection region rej [0, ]2 ,
R= sup
(1 ,2 )rej
|H(1, 2 )|.
(a) Explain why this 2D filter design problem is convex. (b) Find the optimal filter for the specific case with N = 5 and
2 2 rej = {(1 , 2 ) [0, ]2 | 1 + 2 W 2 },
with W = /4. You can approximate R by sampling on a grid of frequency values. Define (p) = (p - 1)/M, p = 1, . . . , M.
(You can use M = 25.) We then replace the exact expression for R above with ^ R = max{|H( (p), (q) )| | p, q = 1, . . . , M, ( (p) , (q) ) rej }. ^ Give the optimal value of R. Plot the optimal frequency response using plot_2D_filt(h), available on the course web site, where h is the matrix containing the coefficients hkl .
2
2. Gini coefficient of inequality. Let x1 , . . . , xn be a set of nonnegative numbers with positive sum, which typically represent the wealth or income of n individuals in some group. The Lorentz curve is a plot of the fraction fi of total wealth held by the i poorest individuals, fi = (1/1T x)
i j=1
x(j) ,
i = 0, . . . , n,
versus i/n, where x(j) denotes the jth smallest of the numbers {x1 , . . . , xn }, and we take f0 = 0. The Lorentz curve starts at (0, 0) and ends at (1, 1). Interpreted as a continuous curve (as, say, n ) the Lorentz curve is convex and increasing, and lies on or below the straight line joining the endpoints. The curve coincides with this straight line, i.e., fi = (i/n), if and only if the wealth is distributed equally, i.e., the xi are all equal. The Gini coefficient is defined as twice the area between the straight line corresponding to uniform wealth distribution and the Lorentz curve: G(x) = (2/n)
n i=1
((i/n) - fi ).
The Gini coefficient is used as a measure of wealth or income inequality: It ranges between 0 (for equal distribution of wealth) and 1 - 1/n (when one individual holds all wealth). (a) Show that G is a quasiconvex function on x Rn \ {0}. + (b) Gini coefficient and marriage. Suppose that individuals i and j get married (i = j) and therefore pool wealth. This means that xi and xj are both replaced with (xi + xj )/2. What can you say about the change in Gini coefficient caused by this marriage?
3
3. Optimal trans-shipment of a commodity. We consider a single commodity, that can be bought or sold in n different locations or markets, at a (given) price pi at location i. Let ui denote the amount sold at location i, where we interpret ui < 0 as meaning that we buy an amount |ui| at location i. The gross revenue from buying and selling the commodity is pT u. At each location i there is a (given) maximum commodity availability ai (which limits how much we can buy there), and a maximum demand for the commodity di (which limits how much we can sell there). Thus, we must have -ai ui di. You can assume that ai and di are nonnegative. (They can be zero, however: di = 0 means that at location i, we can only buy the commodity; we cannot sell it.) We can ship the commodity between the locations, at a cost. Let Sij 0 denote the amount of commodity that we ship from location j to location i. (We can assume that Sii = 0.) The total shipping charge is
n i,j=1
Cij Sij = Tr(C T S),
where Cij are (given) nonnegative shipping rates. Our net profit, including shipping charges, is pT u - Tr(C T S). At location i the total amount shipped out to other locations is
amount received from other locations is
n Sij , so we must have j=1 ui =
n j=1
n
j=1 Sji ,
and the
Sij -
n j=1
Sji.
(In words: the amount of commodity sold at each location is the total amount shipped in to that location, minus the total amount shipped out from the location.) The problem is to choose u and S to maximize the net profit. (a) Suppose there are no shipping charges, i.e., C = 0. Describe a solution, i.e., an optimal u and S. You can assume the prices are in increasing order, i.e., p1 pn . You do not need to justify your answer. Hint. Buy the commodity where it is cheap, and sell it where it is expensive. (b) Solve the instance of the problem given in trans_ship_data.m. Examine the optimal u and S found, and make sure they make sense. (For example: do you buy the commodity at locations where it is cheap, and sell it at locations where it is expensive?) Give a brief informal interpretation of what it means if an entry of the optimal u is zero. Compare the optimal profit obtained to the optimal profit obtained with no shipping costs, as in part (a). (c) Consistent prices. We now remove the limits on buying and selling (i.e., we set all entries in a and d to +). In this case the optimal net profit is either +, 4
or (when 0 the solution is u = 0, S = 0). The first case is called an arbitrage opportunity in economics. In the second case we say the price vector p 0 is consistent, or arbitrage-free. The prices are consistent if no net profit can be made by buying and selling the commodity at different locations and shipping it among the locations. Show that the set of consistent price vectors is a polyhedron. Give the simplest description you can. How would you find the maximum possible consistent price difference between different locations? That is, how you compute the maximum value of |pi - pj | over all consistent price vectors, and all i = j? Find this value for the specific shipping cost data C given in trans_ship_data.m (i.e., ignoring the values of p, a, and d).
5
4. Flux balance analysis in systems biology. Flux balance analysis is based on a very simple model of the reactions going on in a cell, keeping track only of the gross rate of consumption and production of various chemical species within the cell. Based on the known stoichiometry of the reactions, and known upper bounds on some of the reaction rates, we can compute bounds on the other reaction rates, or cell growth, for example. We focus on m metabolites in a cell, labeled M1 , . . . , Mm . There are n reactions going on, labeled R1 , . . . , Rn , with nonnegative reaction rates v1 , . . . , vn . Each reaction has a (known) stoichiometry, which tells us the rate of consumption and production of the metabolites per unit of reaction rate. The stoichiometry data is given by the stoichiometry matrix S Rmn , defined as follows: Sij is the rate of production of Mi due to unit reaction rate vj = 1. Here we consider consumption of a metabolite as negative production; so Sij = -2, for example, means that reaction Rj causes metabolite Mi to be consumed at a rate 2vj . As an example, suppose reaction R1 has the form M1 M2 + 2M3 . The consumption rate of M1 , due to this reaction, is v1 ; the production rate of M2 is v1 ; and the pro duction rate of M3 is 2v1 . (The reaction R1 has no effect on metabolites M4 , . . . , Mm .) This corresponds to a first column of S of the form (-1, 1, 2, 0, . . . , 0). Reactions are also used to model flow of metabolites into and out of the cell. For example, suppose that reaction R2 corresponds to the flow of metabolite M1 into the cell, with v2 giving the flow rate. This corresponds to a second column of S of the form (1, 0, . . . , 0). The last reaction, Rn , corresponds to biomass creation, or cell growth, so the reaction rate vn is the cell growth rate. The last column of S gives the amounts of metabolites used or created per unit of cell growth rate. Since our reactions include metabolites entering or leaving the cell, as well as those converted to biomass within the cell, we have conservation of the metabolites, which can be expressed as Sv = 0. In addition, we are given upper limits on some of the max reaction rates, which we express as v v max , where we set vj = if no upper limit on reaction rate j is known. The goal is to find the maximum possible cell growth rate (i.e., largest possible value of vn ) consistent with the constraints Sv = 0, v 0, v v max .
The questions below pertain to the data found in fba_data.m. (a) Find the maximum possible cell growth rate G , as well as optimal Lagrange multipliers for the reaction rate limits. How sensitive is the maximum growth rate to the various reaction rate limits? (b) Essential genes and synthetic lethals. For simplicity, we'll assume that each re action is controlled by an associated gene, i.e., gene Gi controls reaction Ri . 6
Knocking out a set of genes associated with some reactions has the effect of set ting the reaction rates (or equivalently, the associated v max entries) to zero, which of course reduces the maximum possible growth rate. If the maximum growth rate becomes small enough or zero, it is reasonable to guess that knocking out the set of genes will kill the cell. An essential gene is one that when knocked out reduces the maximum growth rate below a given threshold Gmin . (Note that Gn is always an essential gene.) A synthetic lethal is a pair of non-essential genes that when knocked out reduces the maximum growth rate below the threshold. Find all essential genes and synthetic lethals for the given problem instance, using the threshold Gmin = 0.2G .
7
5. Radiation treatment planning. In radiation treatment, radiation is delivered to a pa tient, with the goal of killing or damaging the cells in a tumor, while carrying out minimal damage to other tissue. The radiation is delivered in beams, each of which has a known pattern; the level of each beam can be adjusted. (In most cases multiple beams are delivered at the same time, in one `shot', with the treatment organized as a sequence of `shots'.) We let bj denote the level of beam j, for j = 1, . . . , n. These must satisfy 0 bj B max , where B max is the maximum possible beam level. The exposure area is divided into m voxels, labeled i = 1, . . . , m. The dose di delivered to voxel i is linear in the beam levels, i.e., di = n Aij bj . Here A Rmn is a (known) matrix j=1 + that characterizes the beam patterns. We now describe a simple radiation treatment planning problem. A (known) subset of the voxels, T {1, . . . , m}, corresponds to the tumor or target region. We require that a minimum radiation dose D target be administered to each tumor voxel, i.e., di D target for i T . For all other voxels, we would like to have di D other , where D other is a desired maximum dose for non-target voxels. This is generally not feasible, so instead we settle for minimizing the penalty E=
iT
((di - D other )+ )2 ,
where ()+ denotes the nonnegative part. We can interpret E as the sum of the squares of the nontarget excess doses. (a) Show that the treatment planning problem is convex. The optimization variable is b Rn ; the problem data are B max , A, T , D target , and D other . (b) Solve the problem instance with data given in the file treatment_planning_data.m. Here we have split the matrix A into Atarget, which contains the rows correspond ing to the target voxels, and Aother, which contains the rows corresponding to other voxels. Give the optimal value. Plot the dose histogram for the target voxels, and also for the other voxels. Make a brief comment on what you see. Re mark. The beam pattern matrix in this problem instance is randomly generated, but similar results would be obtained with realistic data.
8
MIT OpenCourseWare http://ocw.mit.edu
6.079 / 6.975 Introduction to Convex Optimization
Fall 2009
For information about citing these materials or our Terms of Use, visit: http://ocw.mit.edu/terms.