EE364a, Winter 2007-08 Prof. S. Boyd
EE364a Homework 3 solutions
3.42 Approximation width. Let f
0
, . . . , f
n
: R → R be given continuous functions. We
consider the problem of approximating f
0
as a linear combination of f
1
, . . . , f
n
. For
x ∈ R
n
, we say that f = x
1
f
1
+ ··· + x
n
f
n
approximates f
0
with tolerance ǫ > 0 over
the interval [0, T ] if |f(t) − f
0
(t)| ≤ ǫ for 0 ≤ t ≤ T . Now we choose a fixed tolerance
ǫ > 0 and define the approximation width as the largest T such that f approximates
f
0
over the interval [0, T ]:
W (x) = sup{T | |x
1
f
1
(t) + ···+ x
n
f
n
(t) − f
0
(t)| ≤ ǫ for 0 ≤ t ≤ T }.
Show that W is quasiconcave.
Solution. To show that W is quasiconcave we show that the sets {x | W (x) ≥ α} are
convex for all α. We have W (x) ≥ α if and only if
−ǫ ≤ x
1
f
1
(t) + ···+ x
n
f
n
(t) − f
0
(t) ≤ ǫ
for all t ∈ [0, α). Therefore the set {x | W (x) ≥ α } is an intersection of infinitely many
halfspaces (two for each t), hence a convex set.
3.54 Log-concavity of Gaussian cumulative distribution function. The cumulative distribu-
tion function of a Gaussian random variable,
f(x) =
1
√
2π
Z
x
−∞
e
−t
2
/2
dt,
is log-concave. This follows from the general result that the convolution of two log-
concave functions is log-concave. In this problem we guide you through a simple
self-contained proof that f is log-concave. Recall that f is log-concave if and only if
f
′′
(x)f(x) ≤ f
′
(x)
2
for all x.
(a) Verify that f
′′
(x)f(x) ≤ f
′
(x)
2
for x ≥ 0. That leaves us the hard part, which is
to show the inequality for x < 0.
(b) Verify that for any t and x we have t
2
/2 ≥ −x
2
/2 + xt.
(c) Using part (b) show that e
−t
2
/2
≤ e
x
2
/2−xt
. Conclude that
Z
x
−∞
e
−t
2
/2
dt ≤ e
x
2
/2
Z
x
−∞
e
−xt
dt.
(d) Use part (c) to verify that f
′′
(x)f(x) ≤ f
′
(x)
2
for x ≤ 0.
Solution. The derivatives of f are
f
′
(x) = e
−x
2
/2
/
√
2π, f
′′
(x) = −xe
−x
2
/2
/
√
2π.
1
(a) f
′′
(x) ≤ 0 for x ≥ 0.
(b) Since t
2
/2 is convex we have
t
2
/2 ≥ x
2
/2 + x(t −x) = xt − x
2
/2.
This is the general inequality
g(t) ≥ g(x) + g
′
(x)(t − x),
which holds for any differentiable convex function, applied to g(t) = t
2
/2.
Another (easier?) way to establish t
2
/2 ≤ −x
2
/2 + xt is to note that
t
2
/2 + x
2
/2 − xt = (1/2)(x − t)
2
≥ 0.
Now just move x
2
/2 − xt to the other side.
(c) Take exponentials and integrate.
(d) This basic inequality reduces to
−xe
−x
2
/2
Z
x
−∞
e
−t
2
/2
dt ≤ e
−x
2
i.e.,
Z
x
−∞
e
−t
2
/2
dt ≤
e
−x
2
/2
−x
.
This follows from part (c) because
Z
x
−∞
e
−xt
dt =
e
−x
2
−x
.
3.57 Show that the function f(X) = X
−1
is matrix convex on S
n
++
.
Solution. We must show that for arbitrary v ∈ R
n
, the function
g(X) = v
T
X
−1
v.
is convex in X on S
n
++
. This follows from example 3.4.
4.1 Consider the optimization problem
minimize f
0
(x
1
, x
2
)
subject to 2x
1
+ x
2
≥ 1
x
1
+ 3x
2
≥ 1
x
1
≥ 0, x
2
≥ 0.
Make a sketch of the feasible set. For each of the following objective functions, give
the optimal set and the optimal value.
2
(a) f
0
(x
1
, x
2
) = x
1
+ x
2
.
(b) f
0
(x
1
, x
2
) = −x
1
− x
2
.
(c) f
0
(x
1
, x
2
) = x
1
.
(d) f
0
(x
1
, x
2
) = max{x
1
, x
2
}.
(e) f
0
(x
1
, x
2
) = x
2
1
+ 9x
2
2
.
Solution. The feasible set is shown in the figure.
x
1
x
2
(1, 0)
(2/5, 1/5)
(0, 1)
(a) x
⋆
= (2/5, 1/5).
(b) Unbounded below.
(c) X
opt
= {(0, x
2
) | x
2
≥ 1}.
(d) x
⋆
= (1/3, 1/3).
(e) x
⋆
= (1/2, 1/6). This is optimal because it satisfies 2x
1
+x
2
= 7/6 > 1, x
1
+3x
2
=
1, and
∇f
0
(x
⋆
) = (1, 3)
is perpendicular to the line x
1
+ 3x
2
= 1.
4.4 [P. Parrilo] Symmetries and convex optimization. Suppose G = {Q
1
, . . . , Q
k
} ⊆ R
n×n
is
a group, i.e., closed under products and inverse. We say t hat the function f : R
n
→ R
is G-invariant, or symmetric with respect to G, if f (Q
i
x) = f(x) holds for all x and
i = 1, . . . , k. We define
x = (1/k)
P
k
i=1
Q
i
x, which is the average of x over its G-orbit.
We define the fixed subspace of G as
F = {x | Q
i
x = x, i = 1, . . . , k}.
(a) Show that for any x ∈ R
n
, we have
x ∈ F.
(b) Show that if f : R
n
→ R is convex and G- invariant, then f(
x) ≤ f(x).
3
(c) We say the optimization problem
minimize f
0
(x)
subject to f
i
(x) ≤ 0, i = 1, . . . , m
is G-invariant if the objective f
0
is G-invariant, and the feasible set is G-invariant,
which means
f
1
(x) ≤ 0, . . . , f
m
(x) ≤ 0 =⇒ f
1
(Q
i
x) ≤ 0, . . . , f
m
(Q
i
x) ≤ 0,
for i = 1, . . . , k. Show that if the problem is convex and G-invariant, and there
exists an optimal point, then there exists an optimal point in F. In other words,
we can adjoin the equality constraints x ∈ F to the problem, without loss of
generality.
(d) As an example, supp ose f is convex and symmetric, i.e., f(P x) = f(x) for every
permutation P . Show that if f has a minimizer, then it has a minimizer of the
form α1. (This means to minimize f over x ∈ R
n
, we can just as well minimize
f(t1) over t ∈ R.)
Solution.
(a) We first observe that when you multiply each Q
i
by some fixed Q
j
, you get a
permutation of the Q
i
’s:
Q
j
Q
i
= Q
σ(i)
, i = 1, . . . , k,
where σ is a permutation. This is a basic result in group theory, but it’s easy
enough for us to show it. First we note that by closedness, each Q
j
Q
i
is equal
to some Q
s
. Now suppose that Q
j
Q
i
= Q
k
Q
i
= Q
s
. Multiplying by Q
−1
i
on the
right, we see that Q
j
= Q
k
. Thus the mapping from the index i to the index s is
one-to-one, i.e., a permutation.
Now we have
Q
j
x = (1/k)
k
X
i=1
Q
j
Q
i
x = (1/k)
k
X
i=1
Q
σ(i)
x = (1/k)
k
X
i=1
Q
i
x =
x.
This holds for j, so we have
x ∈ F.
(b) Using convexity and invariance of f,
f(
x) ≤ (1/k)
k
X
i=1
f(Q
i
x) = (1/k)
k
X
i=1
f(x) = f(x).
4
(c) Suppose x
⋆
is an optimal solution. Then x
⋆
is feasible, with
f
0
(
x
⋆
) = f
0
((1/k)
k
X
i=1
Q
i
x)
≤ (1/k)
k
X
i=1
f
0
(Q
i
x)
= f
0
(x
⋆
).
Therefore
x
⋆
is also optimal.
(d) Suppose x
⋆
is a minimizer of f. Let
x = (1/n!)
P
P
P x
⋆
, where the sum is over
all permutations. Since
x is invariant under any permutation, we conclude that
x = α1 for some α ∈ R. By Jensen’s inequality we have
f(
x) ≤ (1/n!)
X
P
f(P x
⋆
) = f(x
⋆
),
which shows that
x is also a minimizer.
4.8 Some simple LPs. Give an explicit solution of each of the following LPs.
(a) Minimizing a linear f unction over an affine set.
minimize c
T
x
subject to Ax = b.
Solution. We distinguish three possibilities.
• The problem is infeasible (b 6∈ R(A)). The optimal value is ∞.
• The problem is feasible, and c is orthogonal to the nullspace of A. We can
decompose c as
c = A
T
λ + ˆc, Aˆc = 0.
(ˆc is the component in the nullspace of A; A
T
λ is orthogonal to the nullspace.)
If ˆc = 0, then on the feasible set the objective function reduces to a constant:
c
T
x = λ
T
Ax + ˆc
T
x = λ
T
b.
The optimal value is λ
T
b. All feasible solutions are optimal.
• The problem is feasible, and c is not in the range of A
T
(ˆc 6= 0). The problem
is unbounded (p
⋆
= −∞). To verify this, note that x = x
0
−tˆc is feasible for
all t; as t goes to infinity, the objective value decreases unboundedly.
In summary,
p
⋆
=
+∞ b 6∈ R(A)
λ
T
b c = A
T
λ for some λ
−∞ otherwise.
5
(b) Minimizing a linear function over a halfspace.
minimize c
T
x
subject to a
T
x ≤ b,
where a 6= 0.
Solution. This problem is always feasible. The vector c can be decomposed into
a component parallel to a and a component orthogonal to a:
c = aλ + ˆc,
with a
T
ˆc = 0.
• If λ > 0, the problem is unbounded below. Choose x = −ta, and let t go to
infinity:
c
T
x = −tc
T
a = −tλa
T
a → −∞
and
a
T
x − b = −ta
T
a − b ≤ 0
for large t, so x is feasible for large t. Intuitively, by going very far in the
direction −a, we find feasible p oints with arbitrarily negative objective values.
• If ˆc 6= 0, the problem is unbounded below. Choose x = ba − tˆc and let t go
to infinity.
• If c = aλ for some λ ≤ 0, the optimal value is c
T
ab = λb.
In summary, the optimal value is
p
⋆
=
(
λb c = aλ for some λ ≤ 0
−∞ otherwise.
(c) Minimizing a linear function over a rectangle.
minimize c
T
x
subject to l x u,
where l and u satisfy l u.
Solution. The objective and the constraints are separable: The objective is a
sum of terms c
i
x
i
, each dependent on one variable only; each constraint depends
on only one variable. We can therefore solve the problem by minimizing over
each component of x indepen dently. The optimal x
⋆
i
minimizes c
i
x
i
subject to the
constraint l
i
≤ x
i
≤ u
i
. If c
i
> 0, then x
⋆
i
= l
i
; if c
i
< 0, then x
⋆
i
= u
i
; if c
i
= 0,
then any x
i
in the interval [l
i
, u
i
] is optimal. Therefore, the optimal value of the
problem is
p
⋆
= l
T
c
+
+ u
T
c
−
,
where c
+
i
= max{c
i
, 0} and c
−
i
= max{−c
i
, 0}.
6
(d) Minimizing a linear function over the probability simplex.
minimize c
T
x
subject to 1
T
x = 1, x 0.
What happens if the equality constraint is replaced by an inequality 1
T
x ≤ 1?
We can interpret this LP as a simple portfolio optimization problem. The vector
x represents the allocation of our total budget over different assets, with x
i
the
fraction invested in asset i. The return of each investment is fixed and given by
−c
i
, so our total return (which we want to maximize) is −c
T
x. If we replace the
budget constraint 1
T
x = 1 with an inequality 1
T
x ≤ 1, we have the option of not
investing a portion of the total budget.
Solution. Suppose the components of c are sorted in increasing order with
c
1
= c
2
= ··· = c
k
< c
k+1
≤ ··· ≤ c
n
.
We have
c
T
x ≥ c
1
(1
T
x) = c
min
for all feasible x, with equality if and only if
x
1
+ ··· + x
k
= 1, x
1
≥ 0, . . . , x
k
≥ 0, x
k+1
= ··· = x
n
= 0.
We conclude that the optimal value is p
⋆
= c
1
= c
min
. In the investment interpre-
tation this choice is quite obvious. If the returns are fixed and known, we invest
our total budget in the investment with the highest return.
If we replace the equality with an inequality, the optimal value is equal to
p
⋆
= min{0, c
min
}.
(If c
min
≤ 0, we make the same choice for x as above. Otherwise, we choose
x = 0.)
(e) Minimizing a linear function over a unit box with a total budget constraint.
minimize c
T
x
subject to 1
T
x = α, 0 x 1,
where α is an integer between 0 and n. What happens if α is not an integer (but
satisfies 0 ≤ α ≤ n)? What if we change the equality to an inequality 1
T
x ≤ α?
Solution. We first consider the case of integer α. Suppose
c
1
≤ ··· ≤ c
i−1
< c
i
= ··· = c
α
= ··· = c
k
< c
k+1
≤ ··· ≤ c
n
.
The optimal value is
c
1
+ c
2
+ ··· + c
α
7
i.e., the sum of the smallest α elements of c. x is optimal if and only if
x
1
= ··· = x
i−1
= 1, x
i
+ ··· + x
k
= α − i + 1, x
k+1
= ··· = x
n
= 0.
If α is not an integer, the optimal value is
p
⋆
= c
1
+ c
2
+ ··· + c
⌊α⌋
+ c
1+⌊α⌋
(α − ⌊α⌋).
In the case of an inequality constraint 1
T
x ≤ α, with α an integer between 0 and
n, the optimal value is the sum of the α smallest nonpositive coefficients of c.
4.17 Optimal activity levels. We consider the selection of n nonnegative activity levels,
denoted x
1
, . . . , x
n
. These activities consume m resources, which are limited. Activity
j consumes A
ij
x
j
of resource i, where A
ij
are given. The total resource consumption
is additive, so the total of resource i consumed is c
i
=
P
n
j=1
A
ij
x
j
. (Ordinarily we
have A
ij
≥ 0, i.e., activity j consumes resource i. But we allow the possibility that
A
ij
< 0, which means that activity j actually generates resource i as a by-produ ct.)
Each resource consum ption is limited: we must have c
i
≤ c
max
i
, where c
max
i
are given.
Each activity generates revenue, which is a piecewise-linear concave function of the
activity level:
r
j
(x
j
) =
(
p
j
x
j
0 ≤ x
j
≤ q
j
p
j
q
j
+ p
disc
j
(x
j
− q
j
) x
j
≥ q
j
.
Here p
j
> 0 is the basic price, q
j
> 0 is the quantity discount level, and p
disc
j
is the
quantity discount price, for (the product of) activity j. (We have 0 < p
disc
j
< p
j
.) The
total revenue is the sum of the revenues associated with each activity, i.e.,
P
n
j=1
r
j
(x
j
).
The goal is to choose activity levels that maximize the total revenue while respecting
the resource limits. Show how to formulate this problem as an LP.
Solution. The basic problem can be expressed as
maximize
P
n
j=1
r
j
(x
j
)
subject to x 0
Ax c
max
.
This is a convex optimization problem since the objective is concave and the constraints
are a set of linear inequalities. To transform it to an equivalent LP, we first express
the revenue functions as
r
j
(x
j
) = min{p
j
x
j
, p
j
q
j
+ p
disc
j
(x
j
− q
j
)},
which holds since r
j
is concave. It follows that r
j
(x
j
) ≥ u
j
if and only if
p
j
x
j
≥ u
j
, p
j
q
j
+ p
disc
j
(x
j
− q
j
) ≥ u
j
.
8
We can form an LP as
maximize 1
T
u
subject to x 0
Ax c
max
p
j
x
j
≥ u
j
, p
j
q
j
+ p
disc
j
(x
j
− q
j
) ≥ u
j
, j = 1, . . . , n,
with variables x and u.
To show that this LP is equivalent to the original problem, let us fix x. The last set of
constraints in the LP ensure that u
i
≤ r
i
(x), so we conclude that for every feasible x,
u in the LP, the LP objective is less than or equal to the total revenue. On the other
hand, we can always take u
i
= r
i
(x), in which case the two objectives are equal.
9
Solutions to additional exercises
1. Optimal activity levels. Solve the optimal activity level problem described in exercise
4.17 in Convex Optimization, for the instance with problem data
A =
1 2 0 1
0 0 3 1
0 3 1 1
2 1 2 5
1 0 3 2
, c
max
=
100
100
100
100
100
, p =
3
2
7
6
, p
disc
=
2
1
4
2
, q =
4
10
5
10
.
You can do this by forming the LP you found in your solution of exercise 4.17, or
more directly, using cvx. Give the optimal activity levels, the revenue generated by
each one, and the total revenue generated by the optimal solution. Also, give the
average price per unit for each activity level, i.e., the ratio of the revenue associated
with an activity, to the activity level. (These numbers should be between the basic
and discounted prices for each activity.) Give a very brief story explaining, or at least
commenting on, the solution you find.
Solution. The following Matlab/CVX code solves the p roblem. (Here we write the
problem in a form close to its original statement, and let CVX do the work of refor-
mulating it as an LP!)
A=[ 1 2 0 1;
0 0 3 1;
0 3 1 1;
2 1 2 5;
1 0 3 2];
cmax=[100;100;100;100;100];
p=[3;2;7;6];
pdisc=[2;1;4;2];
q=[4; 10 ;5; 10];
cvx_begin
variable x(4)
maximize( sum(min(p.*x,p.*q+pdisc.*(x-q))) )
subject to
x >= 0;
A*x <= cmax
cvx_end
x
r=min(p.*x,p.*q+pdisc.*(x-q))
10
totr=sum(r)
avgPrice=r./x
The result of the code is
x =
4.0000
22.5000
31.0000
1.5000
r =
12.0000
32.5000
139.0000
9.0000
totr =
192.5000
avgPrice =
3.0000
1.4444
4.4839
6.0000
We notice that the 3rd activity level is the highest and is also the one with the highest
basic price. Since it also has a high discounted price its activity level is higher than the
discount quantity level and it produces the highest contribution to the total revenue.
The 4th activity has a discounted price which is substantially lower then the basic price
and its activity is therefore lower that the discount quantity level. Moreover it require
the use of a lot of resources and therefore its activity level is low.
2. Reformulating constraints in cvx. Each of the following cvx code fragments descr ibes
a convex constraint on the scalar variables x, y, and z, but violates the cvx rule set,
11
and so is invalid. Briefly explain why each fragment is invalid. Then, rewrite each one
in an equivalent form that conforms to the cvx rule set. In your reformulations, you
can use linear equality and inequality constraints, and inequalities constructed using
cvx functions. You can also introduce additional variables, or use LMIs. Be sure to
explain (briefly) why your reformulation is equivalent to the original constraint, if it is
not obvious.
Check your reformulations by creating a small problem that includes these constraints,
and solving it using cvx. Your test problem doesn’t have to be feasible; it’s enough to
verify that cvx processes your constraints without error.
Remark. This looks like a problem about ‘how to use cvx software’, or ‘tricks for
using cvx’. But it really checks whether you understand the various composition rules,
convex analysis, and constraint reformulation rules.
(a) norm( [ x + 2*y , x - y ] ) == 0
(b) square( square( x + y ) ) <= x - y
(c) 1/x + 1/y <= 1; x >= 0; y >= 0
(d) norm([ max( x , 1 ) , max( y , 2 ) ]) <= 3*x + y
(e) x*y >= 1; x >= 0; y >= 0
(f) ( x + y )^2 / sqrt( y ) <= x - y + 5
(g) x^3 + y^3 <= 1; x>=0; y>=0
(h) x+z <= 1+sqrt(x*y-z^2); x>=0; y>=0
Solution.
(a) The lefthand side is correctly identified as convex, b ut equality constraints are only
valid with affine left and right hand sides. Since the norm of a vector is zero if and
only if the vector is zero, we can express the constraint as x+2*y==0; x-y==0, or
simply x==0; y==0.
(b) The problem is that square() can only accept affine arguments, because it is
convex, but not increasing. To correct this use square_pos() instead:
square_pos( square( x + y ) ) <= x - y
We can also reformulate this constraint by introducing an additional variable.
variable t
square( x+y ) <= t
square( t ) <= x - y
Note that, in general, decomposing the objective by introducing new variables
doesn’t need to work. It works in this case because the outer square function is
convex and monotonic over R
+
.
Alternatively, we can rewrite the constraint as
12
( x + y )^4 <= x - y
(c) 1/x isn’t convex, unless you restrict the domain to R
++
. We can write this one
as inv_pos(x)+inv_pos(y)<=1. The inv_pos function has domain R
++
so the
constraints x > 0, y > 0 are (implicitly) included.
(d) The problem is that norm() can only accept affine argument since it is convex
but not increasing. One way to correct this is to introduce new variables u and v:
norm( [ u , v ] ) <= 3*x + y
max( x , 1 ) <= u
max( y , 2 ) <= v
Decomposing the objective by introducing new variables work here because norm
is convex and monotonic over R
2
+
, and in particular over (1, ∞) × (2, ∞).
(e) xy isn’t concave, so this isn’t going to work as stated. But we can express the
constraint as x>=inv_pos(y). (You can switch around x and y here.) Another
solution is to write the constraint as geomean([x,y])>=1. We can also give an
LMI representation:
[ x 1; 1 y ] == semidefinite(2)
(f) This fails when we attempt to d ivide a convex function by a concave one. We can
write this as
quad_over_lin(x+y,sqrt(y)) <= x-y+5
This works because quad_over_lin is monotone descreasing in the second argu-
ment, so it can accept a concave function here, and sqrt is concave.
(g) The function x
3
+ y
3
is convex for x ≥ 0, y ≥ 0. But x
3
isn’t convex for x < 0, so
cvx is going to reject this statement. One way to rewrite this constraint is
quad_pos_over_lin(square(x),x) + quad_pos_over_lin(square(y),y) <= 1
This works because quad_pos_over_lin is convex and increasing in its first ar-
gument, hence accepts a convex function in its first argument. (The function
quad_over_lin, however, is not increasing in its first argument, and so won’t
work.)
Alternatively, and more simply, we can rewrite t he constraint as
pow_pos(x,3) + pow_pos(y,3) <= 1
(h) The problem here is that xy isn’t concave, which causes cvx to reject the state-
ment. To correct this, notice that
q
xy − z
2
=
q
y(x − z
2
/y),
so we can reformulate the constraint as
13
x+z <= 1+geomean([x-quad_over_lin(z,y),y])
This works, since geomean is concave and nondecreasing in each argument. It
therefore accepts a concave function in its first argument.
We can check our reformulations by writing the following feasibility problem in cvx
(which is obviously infeasible)
cvx_begin
variables x y u v z
x == 0;
y == 0;
( x + y )^4 <= x - y;
inv_pos(x) + inv_pos(y) <= 1;
norm( [ u ; v ] ) <= 3*x + y;
max( x , 1 ) <= u;
max( y , 2 ) <= v;
x >= inv_pos(y);
x >= 0;
y >= 0;
quad_over_lin(x + y , sqrt(y)) <= x - y + 5;
pow_pos(x,3) + pow_pos(y,3) <= 1;
x+z <= 1+geomean([x-quad_over_lin(z,y),y])
cvx_end
3. The illumination problem. This exercise concerns the illumination problem described
in lecture 1 (pages 9–11). We’ll take I
des
= 1 and p
max
= 1, so the problem is
minimize f
0
(p) = max
k=1,...,n
|log(a
T
k
p)|
subject to 0 ≤ p
j
≤ 1, j = 1, . . . , m,
(1)
with variable p ∈ R
n
. You will compute several approximate solutions, and compare
the results to the exact solution, for a specific problem instance.
As mentioned in the lecture, the problem is equivalent to
minimize max
k=1,...,n
h(a
T
k
p)
subject to 0 ≤ p
j
≤ 1, j = 1, . . . , m,
(2)
where h(u) = max{u, 1/u} for u > 0. The function h, shown in the figure below, is
nonlinear, nondifferentiable, and convex. To see th e equivalence between (1) and (2),
we note that
f
0
(p) = max
k=1,...,n
|log(a
T
k
p)|
= max
k=1,...,n
max{log(a
T
k
p), log(1/a
T
k
p)}
14
= log max
k=1,...,n
max{a
T
k
p, 1/a
T
k
p}
= log max
k=1,...,n
h(a
T
k
p),
and since the logarithm is a monotonically increasing function, minimizing f
0
is equiv-
alent to minimizing max
k=1,...,n
h(a
T
k
p).
0 1 2 3 4
0
1
2
3
4
5
u
h(u)
The problem instance. The specific problem data are for the geometry shown
below, using the formula
a
kj
= r
−2
kj
max{cos θ
kj
, 0}
from the lecture. There are 10 lamps (m = 10) and 20 patches (n = 20). We take
I
des
= 1 and p
max
= 1. The problem data are given in the file illum_data.m on the
course website. Running this script will construct the matrix A (which has rows a
T
k
),
and plot the lamp/patch geometry as shown below.
0 0.2 0.4 0.6 0.8 1
0
0.5
1
1.5
2
15
Equal lamp powers. Take p
j
= γ for j = 1, . . . , m. Plot f
0
(p) versus γ over
the interval [0, 1]. Graphically determine the optimal value of γ, and the associated
objective value.
You can evaluate the objective function f
0
(p) in Matlab as max(abs(log(A*p))).
Least-squares with saturation. Solve the least-squares problem
minimize
P
n
k=1
(a
T
k
p − 1)
2
= kAp − 1k
2
2
.
If the solution has negative values for some p
i
, set them to zero; if some values are
greater than 1, set them to 1. Give the resulting value of f
0
(p).
Least-squares solutions can be computed using the Matlab backslash operator: A\b
returns the solution of the least-squares problem
minimize kAx − bk
2
2
.
Regularized least-squares. Solve the regularized least-squares problem
minimize
P
n
k=1
(a
T
k
p − 1)
2
+ ρ
P
m
j=1
(p
j
− 0.5)
2
= kAp − 1k
2
2
+ ρkp −(1/2)1k
2
2
,
where ρ > 0 is a parameter. Increase ρ until all coefficients of p are in the interval
[0, 1]. Give the resulting value of f
0
(p).
You can use the backslash operator in Matlab to solve the regularized least-squares
problem.
Chebyshev approximation. Solve the problem
minimize max
k=1,...,n
|a
T
k
p − 1| = kAp − 1k
∞
subject to 0 ≤ p
j
≤ 1, j = 1, . . . , m.
We can think of this problem as obtained by approximating the nonlinear function
h(u) by a piecewise-linear function |u − 1|+ 1. As shown in the figure below, this is a
good approximation around u = 1.
16
0 0.5 1 1.5 2
0
0.5
1
1.5
2
2.5
3
3.5
4
u
|u − 1|+ 1
You can solve the Chebyshev approximation problem using cvx. The (convex) function
kAp −1k
∞
can be expressed in cvx as norm(A*p-ones(n,1),inf). Give the resulting
value of f
0
(p).
Exact solution. Finally, use cvx to solve
minimize max
k=1,...,n
max(a
T
k
p, 1/a
T
k
p)
subject to 0 ≤ p
j
≤ 1, j = 1, . . . , m
exactly. You may find the inv_pos() function useful. Give the resulting (optimal)
value of f
0
(p).
Solution: The following Matlab script finds the approximate solutions using the
heuristic methods proposed, as well as the exact solution.
% illum_sol: finds approximate and exact solutions of
% the illumination problem
clear all;
% load input data
illum_data;
% heuristic method 1: equal lamp powers
% --------------------------------------
nopts=1000;
p = logspace(-3,0,nopts);
f = zeros(size(p));
for k=1:nopts
f(k) = max(abs(log(A*p(k)*ones(m,1))));
17
end;
[val_equal,imin] = min(f);
p_equal = p(imin)*ones(m,1);
% heuristic method 2: least-squares with saturation
% -------------------------------------------------
p_ls_sat = A\ones(n,1);
p_ls_sat = max(p_ls_sat,0); % rounding negative p_i to 0
p_ls_sat = min(p_ls_sat,1); % rounding p_i > 1 to 1
val_ls_sat = max(abs(log(A*p_ls_sat)));
% heuristic method 3: regularized least-squares
% ---------------------------------------------
rhos = linspace(1e-3,1,nopts);
crit = [];
for j=1:nopts
p = [A; sqrt(rhos(j))*eye(m)]\[ones(n,1); sqrt(rhos(j))*0.5*ones(m,1)];
crit = [ crit norm(p-0.5,inf) ];
end
idx = find(crit <= 0.5);
rho = rhos(idx(1)); % smallest rho s.t. p is in [0,1]
p_ls_reg = [A; sqrt(rho)*eye(m)]\[ones(n,1); sqrt(rho)*0.5*ones(m,1)];
val_ls_reg = max(abs(log(A*p_ls_reg)));
% heuristic method 4: chebyshev approximation
% -------------------------------------------
cvx_begin
variable p_cheb(m)
minimize(norm(A*p_cheb-1, inf))
subject to
p_cheb >= 0
p_cheb <= 1
cvx_end
val_cheb = max(abs(log(A*p_cheb)));
% exact solution:
% ---------------
cvx_begin
variable p_exact(m)
minimize(max([A*p_exact; inv_pos(A*p_exact)]))
subject to
p_exact >= 0
18
p_exact <= 1
cvx_end
val_exact = max(abs(log(A*p_exact)));
% Results
% -------
[p_equal p_ls_sat p_ls_reg p_cheb p_exact]
[val_equal val_ls_sat val_ls_reg val_cheb val_exact]
The results are summarized in the following table.
method 1 method 2 method 3 method 4 exact
f
0
(p) 0.4693 0.8628 0.4439 0.4198 0.3575
p
1
0.3448 1 0.5004 1 1
p
2
0.3448 0 0.4777 0.1165 0.2023
p
3
0.3448 1 0.0833 0 0
p
4
0.3448 0 0.0002 0 0
p
5
0.3448 0 0.4561 1 1
p
6
0.3448 1 0.4354 0 0
p
7
0.3448 0 0.4597 1 1
p
8
0.3448 1 0.4307 0.0249 0.1882
p
9
0.3448 0 0.4034 0 0
p
10
0.3448 1 0.4526 1 1
19
