Lecture 11. Matrix Algorithms and Applications Introduction to - - PowerPoint PPT Presentation
Lecture 11. Matrix Algorithms and Applications Introduction to Computer Science, IIIS, Tsinghua Outline } Introduction to divide and conquer } faster matrix multiplication algorithm } Algorithms for matrix operations } LUP decomposition, matrix
} faster matrix multiplication algorithm
} LUP decomposition, matrix inversion, …
} Perfect matching by matrix multiplication } 4-node subgraphs } dominance product / (max,min)-product
} Break up a problem into two sub-problems } Solve each sub-problem recursively } Combine solution to sub-problems to form solution to
} Addition: Given two n-digit integers a and b, compute a + b.
} O(n) bit operations.
} Multiplication. Given two n-digit integers a and b, compute a × b.
} Brute force solution: 𝛪(n2) bit operations.
1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 * 1 1 1 1 1 1 1 + 1 1 1 1 1 1 1 1
} We can divide the n-bit numbers a and b into two parts:
} where a=a1N+a2, b=b1N+b2. (N=2n/2) } so a×b=(a1N+a2)(b1N+b2)=a1b1N2+(a1b2+a2b1)N+a2b2.
} Suppose we can compute the multiplication of two (n/2)-bit
} We need to use it 4 times. } Compute a1b1, a1b2, a2b1, a2b2
} We can divide the n-bit numbers a and b into two parts:
} where a=a1N+a2, b=b1N+b2. (N=2n/2) } so a×b=(a1N+a2)(b1N+b2)=a1b1N2+(a1b2+a2b1)N+a2b2.
} Suppose we can compute the multiplication of two (n/2)-bit
} We need to use it 4 times. } Compute a1b1, a1b2, a2b1, a2b2
T(n) = 4T n/2
recursive calls
1 2 4 3 4 + Θ(n)
add, shift
1 2 3 ⇒ T(n) = Θ(n2)
How can we improve?
} We can divide the n-bit numbers a and b into two parts:
} in fact, a1b1+a1b2+a2b1+a2b2=(a1+a2)(b1+b2). } So if we have already computed a1b1 and a2b2, } a1b2+a2b1=(a1+a2)(b1+b2)-a1b1-a2b2
} We can divide the n-bit numbers a and b into two parts:
} Theorem. [Karatsuba-Ofman, 1962] Can multiply two n-digit
T(n) ≤ T n/2
n/2
recursive calls
1 2 4 4 4 4 4 4 4 3 4 4 4 4 4 4 4 + Θ(n)
add, subtract, shift
1 2 4 3 4 ⇒ T(n) = O(n
log 2 3 ) = O(n1.585 )
T(n) = if n =1 3T(n/2) + n
⎧ ⎨ ⎩
n 3(n/2) 9(n/4) 3k (n / 2k) 3 lg n (2) . . . . . . T(n) T(n/2) T(n/4) T(n/4) T(2) T(2) T(2) T(2) T(2) T(2) T(2) T(2) T(n / 2k) T(n/4) T(n/2) T(n/4) T(n/4) T(n/4) T(n/2) T(n/4) T(n/4) T(n/4)
. . . . . .
T(n) = n 3
2
k k=0 log2 n
∑ =
3 2
1+log2 n −1 3 2 −1
= 3nlog2 3 − 2
} [Karatsuba-Ofman, 1962] O(n1.585)
} log*n is how many times loglog…log n≤1
8 multiplications 4 additions
11 11 11 12 21 12 11 12 12 22 21 21 11 22 21 22 21 12 22 22
1 11 22 11 22 2 21 22 11 3 11 12 22 4 22 21 11 5 11 12 22 6 21 11 11 12 7 12 22 21 22
11 1 4 5 7 12 3 5 21 2 4 22 1 2 3 6
7 multiplications 18 additions/subtractions
11 11 11 12 21 12 11 12 12 22 21 21 11 22 21 22 21 12 22 22
1 11 22 11 22 2 21 22 11 3 11 12 22 4 22 21 11 5 11 12 22 6 21 11 11 12 7 12 22 21 22
11 1 4 5 7 12 3 5 21 2 4 22 1 2 3 6
7 multiplications 18 additions/subtractions
recursive calls
add, subtract
} (using the fast matrix multiplication algorithm)
} (using the fast matrix multiplication algorithm) 可逆矩阵：invertible / nonsingular 不可逆矩阵：not invertible / singular
} L is a unit lower-triangular matrix } U is an upper-triangular matrix } P is a permutation matrix (what’s the inverse of P?)
Permutation matrix: every row and every column has 1 element of 1
⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ ⎡ 1 1 1 1 1 1 1
} For an n×n matrix
} L is a unit lower-triangular matrix } U is an upper-triangular matrix
} Just the Gaussian elimination a wT v A’ 1 v / a
1 1 … 1
a wT
} For an n×n matrix
} L is a unit lower-triangular matrix } U is an upper-triangular matrix
} Just the Gaussian elimination a wT v A’ 1 v / a
1 1 … 1
a wT A’-vwT/a
Recursively perform on A’-vwT/a
} For an n×n matrix
} L is a unit lower-triangular matrix } U is an upper-triangular matrix
} Just the Gaussian elimination a wT v A’ 1 v / a
1 1 … 1
a wT A’-vwT/a
Recursively perform on A’-vwT/a What’s the problem with an arbitrary matrix?
} For an n×n matrix
} L is a unit lower-triangular matrix } U is an upper-triangular matrix
} Just the Gaussian elimination } What if the first element a11=0?
} Change the order of the rows by multiplying a permutation matrix Q
a wT v A’ 1 v / a
1 1 … 1
a wT A’-vwT/a
} Elementary matrix: } Gaussian elimination:
} a sequence of elementary matrix E1, E2, …, Ek s.t. A=E1…EkU } How to deal with the permutation matrix
Left multiply: change rows Right multiply: change columns
a wT v A’ 1 v / a
1 1 … 1
a wT A’-vwT/a
If we change the order of the 2nd line and the i-th line
a wT v A’ 1 v / a
1 1 … 1
a wT A’-vwT/a
If we change the order of the 2nd line and the i-th line
Right multiply: change columns
a wT v A’ 1 v / a
1 1 … 1
a wT A’-vwT/a
If we change the order of the 2nd line and the i-th line
Right multiply: change columns
a wT v A’ 1 v / a
1 1 … 1
a wT A’-vwT/a
If we change the order of the 2nd line and the i-th line
Left multiply: change rows
a wT v A’ 1 v / a
1 1 … 1
a wT A’-vwT/a
If we change the order of the 2nd line and the i-th line
Still a lower triangle matrix!
} LU-decomposition A=LU } When the element aii=0
} Change the order of the rows by multiplying a permutation matrix Q
} Finally A=E1…EkU, for a sequence of elementary matrix } So finally we can write A=PLU
} L is a unit lower-triangular matrix } U is an upper-triangular matrix } P is a permutation matrix
} find inversion } compute determinant
} Solve linear system (Ax=b)
} find inversion } compute determinant } Solve linear system (Ax=b)
} better than O(n3)
} find inversion } compute determinant } Solve linear system (Ax=b)
} So we can recursively compute A-1 and B-1, (also triangle) } Time: T(n)=2T(n/2)+Θ(nω)=Θ(nω)
−1
a wT v A’ 1 v / a
1 1 … 1
a wT A’-vwT/a
} Eliminate elements by order: } O(n3) time.
} Eliminate elements by order:
U X V Y
U X Y
X V Y U
X Y
U
X V Y U
X Y
U
X Y
V
X Y
Recursively run this on Y
X V Y U 1 1 1 1 1 1 1 1
} Inversions:
= −
n i i i I
log 1 1)
} Inversions: } Multiplications:
= −
n i i i I
log 1 1)
= −
n i i i
log 1 1
} Inversions: } Multiplications:
= −
n i i i I
log 1 1)
= −
n i i i
log 1 1
What’s M(n,p) ?
} a trivial method, there is good ways
} Inversions: } Multiplications:
= −
n i i i I
log 1 1)
= −
n i i i
log 1 1
} A reduces to B: we can use an algorithm for B to solve A
} We have shown:
} so LUP-decomposition and matrix inversion has O(nω) time algorithm.
} A reduces to B: we can use an algorithm for B to solve A
} We have shown: } if we can show this reduction, then the three operations has the
} there is an oracle for matrix inversion, compute A×B?
−1
} A reduces to B: we can use an algorithm for B to solve A
} All operation can be solved in O(nω) time.
} A set of vertex-disjoint edges } Matched vertices: the vertices associated with an edge in M } Free vertices: unmatched vertices
} No free vertices
} Problem: decide whether there is a perfect matching in G
} Ai,j=1 iff (i,j)∈E
} Whether there is a perfect matching in G
} Whether there is a perfect matching in G
} σ is the permutation of {1,…,n} } sgn(σ) is +1 or -1
} Whether there is a perfect matching in G
} σ is the permutation of {1,…,n} } sgn(σ) is +1 or -1
} xij is a variable is there is an edge from Ai to Bj } xij=0 if there is no edge from Ai to Bj
} From the definition of determinant, det(X)≡0 if there is no
} The determinant is a polynomial of degree n } We can randomly pick each nonzero xij from {1,…,2n}, and
} if there is no perfect matching, det(X) is always 0 } otherwise, det(X) is not 0 unless we are unlucky enough to pick the
} Computing det(X) takes O(nω) time.
} The determinant is a polynomial of degree n } We can randomly pick each nonzero xij from {1,…,2n}, and
} As in the polynomial identity testing } Computing det(X) takes O(nω) time.
} Proof in the “Theory of Computation” course.
} Transitive Closure: O(nω) } All-pair shortest path (APSP):
} undirected unweighted: O(nω) } directed unweighted: O(n𝜈) } (1+ε)-approximate with weights [0…M]: O(n𝜕 log M/ε)
} Maximum Matching: O(nω) } Dominance Product: O(n(3+ω)/2) } All-pair bottleneck path (APBP):
} vertex-capacitated graphs: O(n𝜈) } edge-capacitated graphs: O(n(3+ω)/2)
} Transitive Closure: O(nω) } All-pair shortest path (APSP):
} undirected unweighted: O(nω) } directed unweighted: O(n𝜈) } (1+ε)-approximate with weights [0…M]: O(n𝜕 log M/ε)
} Maximum Matching: O(nω) } Dominance Product: O(n(3+ω)/2) } All-pair bottleneck path (APBP):
} vertex-capacitated graphs: O(n𝜈) } edge-capacitated graphs: O(n(3+ω)/2)
ω≈2.373 3
(3+ω)/2≈2.687
𝜈≈2.531
} trivial algorithm: checking every triple of vertices (u,v,w) } time: O(n3)
} By FMM?
} (A2)i,j=∑kAi,k·Ak,j
} (A2)i,j=∑kAi,k·Ak,j } (A2)i,j>0 equivalent to ∃k such that Ai,k=Ak,j =1, } that is (i,k), (k,j) ∈E
} (A2)i,j=∑kAi,k·Ak,j } (A2)i,j>0 equivalent to ∃k such that Ai,k=Ak,j =1, } that is (i,k), (k,j) ∈E
1 2 3 4 5 6 7
1 2 3 4 5 6 7
} (A2)2,6=2 sinceA2,7=A7,6 =1, A2,5=A5,6 =1 } Thus time bound is also O(n𝜕)
1 2 3 4 5 6 7
} subgraph on the 4 vertices has 4 edges
1 2 3 4 5 6 7 1 2 3 4 5 6 7
} subgraph on the 4 vertices has 4 edges
1 2 3 4 5 6 7 1 2 3 4 5 6 7
u v
} can only be:
} any edge can be (u,v) } so it is counted 6 times
u v
} with some probability, the number of diamonds will not be a
} we can run this procedure many times, if the graph contains no
2
(u,v)62E
2
(u,v)62E
2
(u,v)62E
} C[i,j]=∑k A[i,k]·B[k,j]
} M[i,j]=(AsB)[i,j]=|{k | A[i,k]≤B[k,j]}|
} M[i,j]=(AsB)[i,j]=|{k | A[i,k]≤B[k,j]}|
} M[i,j]=(AsB)[i,j]=|{k | A[i,k]≤B[k,j]}|
} From V.
2 / ) 1 ( 2 1 − ω
} M[i,j]=(AsB)[i,j]=|{k | A[i,k]≤B[k,j]}|
} M[i,j]=(AsB)[i,j]=|{k | A[i,k]≤B[k,j]}|
} M[i,j]=(AsB)[i,j]=|{k | A[i,k]≤B[k,j]}|
} A[i,1] compare with B[1,j]. } Thus, elements of the i-th column of A are only compared with
} M[i,j]=(AsB)[i,j]=|{k | A[i,k]≤B[k,j]}|
} A[i,1] compare with B[1,j]. } Thus, elements of the i-th column of A are only compared with
} M[i,j]=(AsB)[i,j]=|{k | A[i,k]≤B[k,j]}| } A[i,1] compare with B[1,j]. } Thus, elements of the i-th column of A are only compared
} Thus, define the set Lk={A[i,k]}∪{B[k,j]}, and sort Lk } Partition the sorted list of each Lk into r parts:
} each with 2n/r elements
Increasing
} Ab[i,k]=
} Bb[k,j]=
if i<j, Ai·Bj is always included in the dominance product
} Since the corresponding elements in b-th part of A are ≤ the
} For each b=1…r, compute the Boolean product
} Time: O(r·n𝜕)
} Then we only need to compare elements in the same parts of
} Time: for each part, O(n·(n/r)2), in total: O(rn·(n/r)2)=O(n3/r) } To balance, we let r=n(3-𝜕)/2, total time: O(n(3+𝜕)/2)
} For each b=1…r, compute the Boolean product
} Time: O(r·n𝜕)
} Then we only need to compare elements in the same parts of
} Time: for each part, O(n·(n/r)2), in total: O(rn·(n/r)2)=O(n3/r) } To balance, we let r=n(3-𝜕)/2, total time: O(n(3+𝜕)/2)
} What if A only contains m valid elements?
} other elements of A are +∞
} What if A only contains m valid elements?
} other elements of A are +∞
} For each b=1…r, compute the Boolean product
} Time: O(r·n𝜕)
} Then we only need to compare elements in the same parts of
} Every element in A needs to compare with O(n/r) elements of Bb for
} total time: O(m1/2·n(1+𝜕)/2)
:
p e t s path p ∈ →
} C[i,j]=maxk min{A[i,k], B[k,j]} } which is just the bottleneck path in 3-layer graph:
} G○G ○… ○G = Gn = ((G2)2…)2
n times square about log n times
} C[i,j]=max1≤k≤n{min{A[i,k], B[k,j]}} } So if A[i,k]≤B[k,j], min{A[i,k], B[k,j]}= A[i,k], otherwise it is B[k,j]
} C[i,j]=max1≤k≤n{min{A[i,k], B[k,j]}} } So if A[i,k]≤B[k,j], min{A[i,k], B[k,j]}= A[i,k], otherwise it is B[k,j]
} D[i,j]=max{A[i,k]| A[i,k] ≤ B[k,j]} } D’[i,j]=max{B[k,j]| A[i,k] > B[k,j]} } And, C[i,j]=max{D[i,j], D’[i,j]}
D[i,j]=max{A[i,k]| A[i,k] ≤ B[k,j]}
} Remember we want to compute: } D[i,j]=max{A[i,k]| A[i,k] ≤ B[k,j]} } [Vassilevska,
} Remember we want to compute: } D[i,j]=max{A[i,k]| A[i,k] ≤ B[k,j]} } [Vassilevska,
} Then Ab[i,j]=
} Then there exists a Ai,k in Ap such that Ai,k≤Bk,j, so Di,j is in Ap.
} Then there exists a Ai,k in Ap-1 such that Ai,k≤Bk,j, so Di,j is in Ap-
1.
} Remember we want to compute: } D[i,j]=max{A[i,k]| A[i,k] ≤ B[k,j]} } [Vassilevska,
} Then Ab[i,j]=
} Compute the dominance product AbsB for b=1…s } Then for each pair of (i,j), find the maximum b such that
} Remember we want to compute: } D[i,j]=max{A[i,k]| A[i,k] ≤ B[k,j]} } [Vassilevska,
} Compute the dominance product AbsB for b=1…s
} Time: O(m1/2·n(1+𝜕)/2) for each sparse dominance product where
} Then for each pair of (i,j), find the maximum b such that
} Time: O(n/s) for each i,j
} T
ω≈2.376 3
(3+ω)/2 2+ω/3
Algorithms for bottleneck paths Running Time Source Edge- weighted Single-source O(m+nlogn) Single-source, single- destination O(m log*n) Gabow & Tarjan All-pair (trivial) O(mn) All-pair (2007) O(n2+ω/3)≈O(n2.792) Vassilevska,William s and Yuster All-pair (2009) O(n(3+ω)/2)≈O(n2.688) Duan, Pettie Vertex- weighted All-pair(2007) O(nμ)≈O(n2.575) Shapira, Yuster and Zwick
} Algorithms for matrix operations of bound O(n2.373):
} matrix multiplication } inversion } LUP-decomposition } solving linear equations } determinant } but not eigenvalues!
} Applications:
} Perfect matching in bipartite graphs } 4-node subgraphs } dominance product / (max,min)-product