Algorithm Design and Analysis
LECTURE 16
Dynamic Programming
(plus FFT Recap)
Adam Smith
9/24/2008
A. Smith; based on slides by E. Demaine, C. Leiserson, S. Raskhodnikova, K. Wayne
Discrete Fourier Transform
Coefficient to point-value. Given a polynomial a0 + a1 x + ... + an-1 xn-1,
evaluate it at n distinct points x0, ... , xn-1.
Key idea: choose xk = k where is principal nth root of unity.
2 = 1 = i
3
4 = 2 = -1
1
n=8
7
5
6 = 3 = -i
Discrete Fourier transform
Fourier matrix Fn
0 = 0 = 1
Fast Fourier Transform
Goal. Evaluate a degree n-1 polynomial A(x) = a0 + ... + an-1 xn-1 at its nth
roots of unity: 0, 1, , n-1.
Divide. Break polynomial up into even and odd powers.
Aeven(x) = a0 + a2x + a4x2 + + an/2-2 x(n-1)/2.
Aodd (x) = a1 + a3x + a5x2 + + an/2-1 x(n-1)/2.
A(x) = Aeven(x2) + x Aodd(x2).
Conquer. Evaluate degree Aeven(x) and Aodd(x) at
set of all squares of nth roots of unity
= set of all nth roots of unity: 0, 1, , n/2-1
Combine.
A(k+n) = Aeven(k) + k Aodd(k), 0 k < n/2
A(k+n/2) = Aeven(k) - k Aodd(k), 0 k < n/2
k+n/2 = -k
k = (k)2 = (k+n/2)2
When n is even,
there are only
n/2 possible
values of 2k
FFT Algorithm
fft(n, a0,a1,,an-1) {
if (n == 1) return a0
(e0,e1,,en/2-1) FFT(n/2, a0,a2,a4,,an-2)
(d0,d1,,dn/2-1) FFT(n/2, a1,a3,a5,,an-1)
for k =
k
yk+n/2
yk+n/2
}
0 to n/2 - 1 {
e2ik/n
e k + k d k
e k - k d k
return (y0,y1,,yn-1)
}
FFT Summary
Theorem. FFT algorithm evaluates a degree n-1 polynomial at each of
the nth roots of unity in O(n log n) steps.
assumes n is a power of 2
Running time. T(2n) = 2T(n) + O(n) T(n) = O(n log n).
O(n log n)
coefficient
representation
point-value
representation
Point-Value to Coefficient Representation: Inverse DFT
Goal. Given the values y0, ... , yn-1 of a degree n-1 polynomial at the n
points 0, 1, , n-1, find unique polynomial a0 + a1 x + ... + an-1 xn-1 that
has given values at given points.
Inverse DFT
Fourier matrix inverse (Fn)-1
Inverse FFT
Claim. Inverse of Fourier matrix is given by following formula.
Consequence. To compute inverse FFT, apply same algorithm but use
-1 = e -2 i / n as principal nth root of unity (and divide by n).
Polynomial Multiplication
Theorem. Can multiply two degree n-1 polynomials in O(n log n) steps.
coefficient
representation
FFT
coefficient
representation
O(n log n)
inverse FFT
point-value multiplication
O(n)
O(n log n)
Integer Multiplication
Integer multiplication. Given two n bit integers a = an-1 a1a0 and
b = bn-1 b1b0, compute their product c = a b.
Convolution algorithm.
Form two polynomials.
Note: a = A(2), b = B(2).
Compute C(x) = A(x) B(x).
Evaluate C(2) = a b.
Running time: O(n log n) complex arithmetic steps.
Theory. [Schnhage-Strassen 1971] O(n log n log log n) bit operations.
[Martin Frer (Penn State) 2007] O(n log n 2log* n) bit operations.
Practice. [GNU Multiple Precision Arithmetic Library] GMP proclaims
to be "the fastest bignum library on the planet." It uses brute force,
Karatsuba, and FFT, depending on the size of n.
Wrap-up Divide and Conquer
Look for recursive structure
Steps:
Divide: Break into smaller subproblems
Conquer: Solve subproblems
Combine: compute final solution
Often helps to compute extra information in
subproblems
e.g. # inversions: recursive call also sorts each piece
If pieces independent, get parallelization
Chapter 5
Dynamic Programming
Design Techniques So Far
Greedy. Build up a solution incrementally,
myopically optimizing some local criterion.
Divide-and-conquer. Break up a problem into
subproblems, solve subproblems, and combine
solutions.
Dynamic programming. Break problem into
overlapping subproblems, and build up solutions
to larger and larger sub-problems.
10/8/2007
S. Raskhodnikova; based on slides by K. Wayne.
Dynamic Programming History
Bellman. [1950s] Pioneered the systematic study of dynamic programming.
Etymology.
Dynamic programming = planning over time.
Secretary of Defense was hostile to mathematical research.
Bellman sought an impressive name to avoid confrontation.
"it's impossible to use dynamic in a pejorative sense"
"something not even a Congressman could object to"
Reference: Bellman, R. E. Eye of the Hurricane, An Autobiography.
Dynamic Programming Applications
Areas.
Bioinformatics.
Control theory.
Information theory.
Operations research.
Computer science: theory, graphics, AI, compilers, systems, .
Some famous dynamic programming algorithms.
Unix diff for comparing two files.
Viterbi for hidden Markov models / decoding convolutional codes
Smith-Waterman for genetic sequence alignment.
Bellman-Ford for shortest path routing in networks.
Cocke-Kasami-Younger for parsing context free grammars.
Fibonacci Sequence
Sequence defined by
a1 = 1
a2 = 1
an = an-1 + an-2
1, 1, 3, 5, 8, 13, 21, 34,
How should you compute the Fibonacci sequence?
Recursive algorithm:
Running Time?
Fib(n)
1. If n =1 or n=2, then
2.
return 1
3. Else
4.
a = Fib(n-1)
5.
b = Fib(n-2)
6.
return a+b
Review Question
Prove that the solution to the recurrence T(n)=T(n-1)+T(n-2) + (1) is
exponential in n.
Computing Fibonacci Sequence Faster
Observation: Lots of redundancy! The recursive algorithm only solves
n-1 different sub-problems
Memoization: Store the values returned by recursive calls in a subtable
Resulting Algorithm:
Fib(n)
1. If n =1 or n=2, then
2.
return 1
3. Else
4.
f[1]=1; f[2]=1
5.
For i=3 to n
6.
f[i] f[i-1]+f[i-2]
7.
return a+b
Linear time!
(In fact, a log(n)-time algorithm exists.)
Weighted Interval Scheduling
Weighted interval scheduling problem.
Job j starts at sj, finishes at fj, and has weight or value vj .
Two jobs compatible if they don't overlap.
Goal: find maximum weight subset of mutually compatible jobs.
a
b
c
d
e
f
g
h
0
Time
10
Unweighted Interval Scheduling Review
Recall. Greedy algorithm works if all weights are 1.
Consider jobs in ascending order of finish time.
Add job to subset if it is compatible with previously chosen jobs.
Observation. Greedy algorithm can fail spectacularly if arbitrary
weights are allowed.
weight = 999
weight = 1
a
Time
0
10
11
Weighted Interval Scheduling
Notation. Label jobs by finishing time: f1 f2 . . . fn .
Def. p(j) = largest index i < j such that job i is compatible with j.
Ex: p(8) = 5, p(7) = 3, p(2) = 0.
1
2
3
4
5
6
7
8
0
Time
10
11
Dynamic Programming: Binary Choice
Notation. OPT(j) = value of optimal solution to the problem consisting
of job requests 1, 2, ..., j.
Case 1: OPT selects job j.
collect profit vj
can't use incompatible jobs { p(j) + 1, p(j) + 2, ..., j - 1 }
must include optimal solution to problem consisting of remaining
compatible jobs 1, 2, ..., p(j)
optimal substructure
Case 2: OPT does not select job j.
must include optimal solution to problem consisting of remaining
compatible jobs 1, 2, ..., j-1
Weighted Interval Scheduling: Brute Force
Brute force algorithm.
Input: n, s1,,sn
f1,,fn
v1,,vn
Sort jobs by finish times so that f1 f2 ... fn.
Compute p(1), p(2), , p(n)
Compute-Opt(j) {
if (j = 0)
return 0
else
return max(vj + Compute-Opt(p(j)), Compute-Opt(j-1))
}
Weighted Interval Scheduling: Brute Force
Observation. Recursive algorithm fails spectacularly because of
redundant sub-problems exponential algorithms.
Ex. Number of recursive calls for family of "layered" instances grows
like Fibonacci sequence.
5
2
3
3
4
5
p(1) = 0, p(j) = j-2
Weighted Interval Scheduling: Memoization
Memoization. Store results of each sub-problem in a table;
lookup as needed.
Input: n, s1,,sn
f1,,fn
v1,,vn
Sort jobs by finish times so that f1 f2 ... fn.
Compute p(1), p(2), , p(n)
for j = 1 to n
M[j] = empty
M[0] = 0
global array
M-Compute-Opt(j) {
if (M[j] is empty)
M[j] = max(wj + M-Compute-Opt(p(j)), M-Compute-Opt(j-1))
return M[j]
}
Weighted Interval Scheduling: Running Time
Claim. Memoized version of algorithm takes O(n log n) time.
Sort by finish time: O(n log n).
Computing p() : O(n log n) via sorting by start time.
M-Compute-Opt(j): each invocation takes O(1) time and either
(i) returns an existing value M[j]
(ii) fills in one new entry M[j] and makes two recursive calls
Progress measure = # nonempty entries of M[].
initially = 0, throughout n.
(ii) increases by 1 at most 2n recursive calls.
Overall running time of M-Compute-Opt(n) is O(n).
Remark. O(n) if jobs are pre-sorted by start and finish times.
Weighted Interval Scheduling: Finding a Solution
Q. Dynamic programming algorithms computes optimal value.
What if we want the solution itself?
A. Do some post-processing.
Run M-Compute-Opt(n)
Run Find-Solution(n)
Find-Solution(j) {
if (j = 0)
output nothing
else if (vj + M[p(j)] > M[j-1])
print j
Find-Solution(p(j))
else
Find-Solution(j-1)
}
# of recursive calls n O(n).
Weighted Interval Scheduling: Bottom-Up
Bottom-up dynamic programming. Unwind recursion.
Input: n, s1,,sn
f1,,fn
v1,,vn
Sort jobs by finish times so that f1 f2 ... fn.
Compute p(1), p(2), , p(n)
Iterative-Compute-Opt {
M[0] = 0
for j = 1 to n
M[j] = max(vj + M[p(j)], M[j-1])
}
