7High-Performance General Solver
for Extremely Large-Scale Semidefinite
Programming Problems
The multiple GPU calculation
for the classical spin model
Knowledge Acquisition froma Large Web Corpus and its Applications

High-Performance General Solver for Extremely
Large-Scale Semidefinite Programming Problems
Katsuki Fujisawa* Toshio Endo***Chuo University & JST CREST**Tokyo Institute of Technology & JST CREST

Semidef inite programming (SDP) is one of the most impor tant problems among optimiz ation problems atpresent. It is relevant to a wide range of f ields such as combinatorial optimization, struc tural optimization,
control theory, economics, quantum chemistry, sensor network location and data mining. The capability to solve
extremely large-scale SDP problems will have a significant effect on the current and future applications of SDP. In
1995, Fujisawa et al. started the SDPA (Semidefinite programming algorithm) Project aimed at solving large-scale
SDP problems with high numerical stability and accuracy. SDPA is one of the main codes to solve general SDPs.
SDPARA is a parallel version of SDPA on multiple processors with distributed memory, and it replaces two majorbottleneck parts (the generation of the Schur complement matrix and its Cholesky factorization) of SDPA by theirparallel implementation. In particular, it has been successfully applied to combinatorial optimization and truss
topology optimization. The new version of SDPARA (7.5.0-G) on a large-scale supercomputer called TSUBAME 2.0 [12]at the Tokyo Institute of Technology has successfully been used to solve the largest SDP problem (which has over 1.48
million constraints), and created a new world record. Our implementation has also achieved 533 TFlops in double
precision for large-scale Cholesky factorization using 2,720 CPUs and 4,080 GPUs. [1]Keywords : GPGPU, Dense Matrix Algebra, Semidefinite Programming, General Optimization Solver
matrices: the computation of the Schur complement matrix

1

Introduction(SCM) and the Cholesky factorization of the SCM.
The semidefinite programming algorithm parallelversion (SDPARA) [5], [6] is a parallel implementation of the SDPA
library for multiple processor systems with distributed memory.

Semidefinite Programming (SDP) is a subfield of mathematicalIn this implementation, performance bottlenecks are relaxed

programming and is used to optimize an objective function overby parallelization techniques using MPI, OpenMP, ScaLAPACK

positive semidefinite matrices. SDP has been regarded as one of[7]

the most important optimization problems for several reasons.

the input data and stores them and all variables in the process

Solving extremely large-scale SDP problems is considered tomemor y space, while the SCM data are divide d b et we en

be important and challenging. For example, it would be usefulprocesses. Thus, SDPARA can compute each row of the SCM into obtain the optimal solution of the quadratic assignment

parallel, and apply the parallel Cholesky factorization provided

problem (QAP), one of the most difficult NP-hard problems. If weby ScaLAPACK to the SCM. We have repor ted that SDPARA6, BLACS, optimized BLAS, and MUMPS [8]. Each process reads

could solve an SDP problem with more than 10 constraints, weachieves high scalability for large-scale SDP problems [2], [5], [6]could obtain a strong lower bound of the QAP with a problem

on supercomputers and have confirmed that SDPARA is muchsize greater than or equal to 40.

faster than other parallel implementations such as PCSDP [9] andFor solving SDP problems in polynomial time, theprimal-dual interior-point method (PDIPM)[11]PDSDP [10] when solving large-scale sparse SDP problems [5], [6].

is well known as

By using existing SDPARA (version 7.4.0), we have successfully

a powerful and stable method. In 1995, Fujisawa et al. started

solved large-scale SDP problems with a sparse SCM and a size of

the SDPA Project[4]

aimed at solving large-scale SDP problems

about 300,000, and created a world record in 2011 [2].

with high numerical s tabilit y and accurac y. Semidef inite
programming algorithm (SDPA)

02

[4]

T he target of this study is highly challenging: inis one of the main libraries

order to accommodate problems related to the QAP and truss

f o r s o l v i n g g e n e r a l SD P p r o b l e m s t h a t h av e n o s p e c i a l

topology optimization, our objective is to solve SDP problems

struc tures such as a rank-1 condition, and an optimization

with a dense SCM and a size larger than 10 6 . Aside from an

sof tware package is used along with the PDIPM for solving

increase in problem sizes, changes in the characteristics ofthe standard SDP problem (formulation (1) in Section 2). SDPA

problems (the density of the SCM) also pose new challenges. Let

incorporates special data structures for handling block diagonal

m be the SCM size (which equals the number of constraints) anddata matrices, and efficient techniques for computing searchn be the size of the data matrix. For the first bottleneck part

directions when problems become large scale and/or sparse(computation of the SCM), we have previously developed an[3]

algorithm that decreases its time complexity O(mn +m n ) to[4]

. SDPA achieves good performance on a single processor ;

3

2

[5]

2

2

however, when the size of an SDP problem increases, the use

O(m ) by exploiting the sparsity of the data matrixof the PDIPM introduces two major performance bottlenecks

if the SCM is large and dense, the second bot tleneck par t. However,in each iteration, even if we exploit the sparsity of the data

(Cholesky factorization of the SCM), whose complexity is O(m ) ,3

dominates the entire execution time. In this study, we accelerate

Note that

this part by using massively parallel GPUs with computational

share the same block structure. We define n max as max｛ n 1 ,..., n h ｝.

performance much higher than that of CPUs. In order to achieve

For the block s wherescalable performance with thousands of GPUs, we utilize a high-semidefiniteness are equivalent to the cons traint s of theperformance BLAS kernel along with optimization techniques

nonnegative or thant. Such block s are sometimes called LP

to overlap computation, PCI-Express communication, and MPI

blocks.

communication.

The size of given SDP problems can be roughly measured inThis paper describes the design and implementationof a new version of SDPARA (version 7.5.0-G). We have selected
three types of large-scale SDP instances associated with recent

and the variable matrices

and, the constraints of positive

terms of five numbers:(1) m : the number of equality constraints in the dual form
(which equals the size of the SCM)

applications of SDP in truss topology design, combinatorial(2) n : the size of the variable matrices

optimiz ation, and quantum chemis tr y; however, SDPAR A(3) n max : the size of the largest block of input data matrices

can also achieve high performance when solving large-scale(4) nnz : the total number of nonzero elements in all data

SDP problems in other important application areas. We haveandmatrices

conduc ted the per formance evaluation of SDPARA 7.5.0 - G

Figure 1 shows the largest SDP problem, which is solved

o n T SUBA ME 2. 0, w hi ch is a p e t a - s c al e GPU -a cce l er ate din Section 6, and its block diagonal structure. The parameters m , n ,supercomputer at the Tok yo Institute of Technology. In the

n max , and nnz in this SDP problem are 1,484,406, 1,779,204, 1,682,

evaluation, we solved the largest SDP problem (which has overand 23,476,376, respectively.

1.48 million constraints), and created a new world record in 2012.
Our implementation also achieved 533 TFlops in double precision
for large-scale Cholesky factorization using 4,080 GPUs.

2

Semidefinete programming

The standard form SDP is the following primal-dual form.

Figure 1 Largest SDP problem and its block
diagonal structure
We usen

for the space of n×n symmetric matrices. The notation
indicates thatnis a positive semidefinite(positive definite) matrix. The innerproduct between
and

n

n

is defined by

Basic framework of the primal-dual
interior-point method

3

In most SDP applications, it is common for the input
data matrices

0

,...,

m

to share the same diagonal blo ck

structure ( n1 ,..., n h ). Each input data matrix

K( k = 1 ,...,m ) consistsof sub-matrices in the diagonal positions, as follows:

In this section, we explain the basic framework of the PDIPM
on which SDPA and SDPARA are based. The most important
theoretical aspects of the PDIPM is that it solves both primal anddual forms simultaneously and finds their optimal solution inpolynomial time.
The Karush-Kuhn-Tucker (KKT) conditions theoreticallyguarantee that a point (*,*,*) satisfying the system in (2)below is an optimal solution of (1) when the so-called Slater’s
condition is satisfied.

03

High-Performance General Solver for Extremely
Large-Scale Semidefinite Programming Problems

As we have shown in Figure 2, in the PDIPM, the algorithm

Step 3 : Maximize the step sizes ( lengths ) such that thefollowing positive definiteness conditions are satisfied.

starts from a feasible or infeasible point. In each iteration, it
computes the search direction ( d , d,d) from the current

point towards the optimal solution, decides the step size, andadvances by the step size in the search direction. If the current
point reaches a small neighborhood of the optimal solution,

Step 4 : Update the current point as follows:

the PDIPM terminates the iteration and returns the approximate
optimal solution. The PDIPM is described in many papers (See [11] ).
The framework we use relies on the HRVW/KSH/M approach [11], andSet

we will use the appropriate norms ||・|| for matrices and vectors.

, and return to Step 1.

Steps 2a and 2b correspond to the first and second
bottleneck parts defined in Section 1, respectively. We shall define

Algorithmic Framework of PDIPM

the computation in Steps 2a and 2b as ELEMENTS and CHOLESKY,
respectively. As indicated in [5] , ELEMENTS (O (mn + m n ) ) andStep 0 : Choose a feasible or infeasible initial point3such that, the boundary, the threshold parameter,and the iteration number s = 0 .

total execution time of the PDIPM. Therefore, researchers havefocused on reducing the time taken for these steps [4]. We have3

2

2

reduced the time complexity of ELEMENTS from (O (mn + m n ) ) to

Step 1 : Compute the residuals of the primal feasibility
dual feasibility

2

CHOLESKY (O (m ) ) have often accounted for 80% to 90% of the

Set the centering parameterparameter

2

3

, the

O ( m ) when solving the sparse SDP problem [3].
2, and the primal-dual gap g :If( namely, all residuals above

are sufficiently small ), stop the iteration and output
(

s

,

s

,

s

) as an approximate optimal solution.

Step 2 : Compute the the search direction ( d , dStep 2a : Compute the SCM,d).

by the formula

Step 2b : Apply the Cholesky factorization toa lower triangular matrix

and obtainsuch that

Step 2c : Obtain a component of the search direction d.
by
Figure 2 Primal-dual interior-point method

solving the equations
for the right-hand-side vector r computed as

Accelerating cholesky Factorization
where

with

Step 2d : Compute the remaining components of the search
direction ( d4,d) as follows:

A. Challenges

As described above, the total execution timeof SDPAR A is dominate d by that of CHOLE SK Y (Cho l esk yfactorization of SCM

04) whenis sufficiently large and dense.

Generally, dense matrix computation can be significantly

kernel dominates the execution time. Note that updating

accelerated by harnessing GPGPU computing; for example,

the lower-triangular part is sufficient; thus, we can omit theEndo et al. have shown that the performance of the Linpack

computation of the unused upper part.
T he b asic appro aches we app lie d to acceler ate thisbenchmark (LU factorization with pivoting) can be scalably
increased by using more than 1200 accelerators [13]. The keys toscalability include

algorithm are as follows.
1) We invoke one MPI process per GPU to drive it, and thus,

・ Overlapping computation on GPUs and PCI-Express (PCIe)three processes per node are invoked on TSUBAME 2.0 nodes.
2) It is important to use the fast DTRSM and DGEMM BLAS

communication between GPUs and CPUs.
・ Configuring the block size nb so that data reuse is promoted

kernels on GPUs. We used highly tuned kernels developed

and the amount of PCIe communication is reduced. In thisby NVIDIA, which are faster than those in the off icial CUBLAS

paper, nb = 1024 .

library. The DGEMM kernel achieves about 350 GFlops (on-While the parallel Cholesky factorization algorithm andboard speed) on a M2050 GPU, while CUBLAS DGEMM
achieves about 300 GFlops 1.

Linpack have many common points, the former poses new3) On GPU clusters, we have to decide where the datastructure

challenges because of the following differences.
・ In each computation step, the part of the matrix to be

is located since the GPU device memory is separated fromthe host memory. In order to accommodate larger sizes of

updated is the lower-triangle part, rather than a rectangle.

, we store it on the host memory.

Upon two-dimensional block-cyclic distribution, the shapeof the updated part in each process becomes more complex.

Approaches 2) and 3) indicate that we need to divide the・ (Related to the abaove difference) The computation amountmatrices into parts smaller than the device memory capacity

per step is halved compared to LU factorization. This makes

and send the input matrices to the GPU via PCIe in order tothe computation / communication ratio even worse.

perform partial computation. Of course, without overlapping

・ Cholesky factorization requires additional work called

computation and communication (as in “version 1” in Figure 3),”panel transposition.”the performance is strictly restricted. To ensure high performance,
the following optimization and configuration are adopted.
・ When the size of the partial matrix to be updated by a single

B. Implementation and Optimization

Hereafter, we call SDPARA whose CHOLESKY component

GPU is r × s , the computation cost in the “update” phase isis a cce l e r a te d b y G PUs “ SD PA R A (v e r s i o n 7. 5. 0 - G ).” O u r

O (r・s・nb) , while the communication cost is O(r・s +r・nb+s・nb) .

accelerated CHOLESKY component has properties similar to

To reduce the relative communication cost, the block size

the pdpotrf function of ScaLAPACK. The dense matrix, withnb should be sufficiently large. Here there is a trade-off sincea size of m × m , is distributed among MPI processes in the two-very large nb degrades load balancing. After a preliminary
evaluation, we set nb to be 1024.

dimensional block-cyclic distribution with block size n b . When
we let mb =⎡m /nb ⎤be the number of blocks that are aligned

・ In order to reduce and hide the PCIe communication cost,in a row or a column, the CHOLESKY algorithm consists of mb

we overlap GPU computation and PCIe communication.
With these methods (“version 2” in the figure) we can enjoy

steps. A single (k -th) step proceeds as follows:・ Diagonal block factorization: The k - th diagonal block is
Cholesky-factorized locally. Then, the result block isinter-node MPI communication costs still restrict the scalability.

broadcast to processes in the k - th process column.

Although each node in TSUBAME 2.0 has a wide injec tion

・ Panel factorization: The k - th block columns are calledbandwidth (8 GB/s), the communication costs increase relatively“panel”, and the panel is factorized by using the dtrsm

・ Panel broadcast and transposition: We need to broadcast
row-wise, obtain the transposition oftwhen the computation is accelerated. Here, we further promote theoverlapping policy described above; we overlap all computations,

BLAS kernel.
, and broadcast

PCIe communication, and MPI communication (“version 3” inthe figure). For this, the “Panel broadcast and transposition”
and “Update” phases are reorganized; the transposed panel

column-wise.

t・ Update: This is the most computation-intensive part. Eachis divided into pieces before broadcasting, and each processprocess updates its own part of the rest matrix, taking thetransfers them to the GPU just after partial broadcast is completed.

corresponding part ofThen
1

the accelerated performance. However, we have noticed that

'=

'-

×

t

and

t

. Let

As future improvement, we could overlap MPI communication' be the rest matrix.

is computed. Thus, the DGEMM BLAS

forand computation by using an optimization technique

Unfortunately, the speed of 350 GFlops is still far from the theoretical peak of 515 GFlops. NVIDIA staff mention that this is because of the limitation of the architectureof Fermi generation GPUs and that the speed will increase in next-generation GPUs.

05

High-Performance General Solver for Extremely
Large-Scale Semidefinite Programming Problems

c all e d “l o o k ahea d,” w hich has b e en intro duce d in High Performance Linpack [14].

of nodes; the difference is 10% for 64 nodes.
Figure 5 shows the performance of CHOLESKY on almost
all nodes of TSUBAME 2.0 (1,360 nodes with 4,080 GPUs). For the

Performance evaluation

5

QAP8 problem, where the SCM size m is 1,484,406, we achieve

533TFlops, which is the maximum speed achieved to date for amathematical optimization problem. As in the previous graph,
we see that a larger problem corresponds to a higher speed.
Since the peak performance of 1,360 nodes is 2,293 TFlops,

We have conducted the performance evaluation of SDPARA

the relative performance with QAP8 is 24%. This is lower than

7.5.0-G with GPU-accelerated Cholesk y factorization on the

results shown in the previous graph. However, we note that there

TSUBAME 2.0 supercomputer. Table I shows the problems

is plenty of space in the node memory owing to the optimizationconsidered in the evaluation; the problems have characteristics

for memory usage reduction. In the case of QAP8 on 1,360 nodes,that lead to the SCM being dense. We used NVIDIA CUDA 4.0,

only 21 GB of the 54 GB memory is used per node. Thus, we can

MVAPICH2 1.5.1, and Intel compiler 11.1 as the underlying system

accommodate larger a problem with m 〜 2, 000, 000 , whichsoftware. We used BLAS kernels from NVIDIA for use on GPUs, as

would further increase the speed.

mentioned in SectionIV-B.
We use different MPI processes for the different GPUs;
thus, three processes run on a TSUBAME 2.0 node. A single
process is configured to drive a GPU and three CPU cores, byconsidering the CPU affinity by using cudaSetDevice CUDA API
and sched_setaffinity systemcall; for example, a process thatuses GPU 0 in each node is bound to the cores in CPU 0. Three
CPU cores are used in INIT, ELEMENTS, and so on, while they
are unused in the DGEMM or DTRSM kernels in CHOLESKY. For

Version 1 No overlappingthese kernels, only GPUs are used 2 .
A. Performance of Accelerated Cholesky Factorization

First, the measured performance for a moderate numberof nodes is shown in Figure 4. The graph shows the speed of the
CHOLESKY component in TeraFlops. We see that its speed for
the QAP5 problem reaches 46.5 TFlops for 64 nodes (192 GPUs).
This high performance is obviously impossible without GPUs
since the theoretical peak performance of CPUs on 64 nodes is
64 × 0.1408 = 9.01 TFlops. Compared to the peak performanceof CPUs and GPUs (6 4 ×1.686 = 107.9 Tf lops), the relative

Version 2 GPU computation and PCIe
communication are overlapped

performance is 43.1%. This gap results from several factors, andthe most important one is the limited performance of on-board
DGEMM (only 68% of peak), which will be improved in nextgeneration Kepler GPUs. A comparison between QAP5 and QAP4shows that the former shows better Flop performance sincea larger matrix size is advantageous in CHOLESKY with O(m 3)computation and O(m 2) communication.
T he graph also shows the ef fec t of overlapping MPI
communication by presenting a comparison bet ween t wo
CHOLESKY implementations: Version 3 and Version 2 describedin Section 4 - B. By comparing QAP4/v3 and QAP4/v2 , we can seethat the effect of overlapping increases with the larger number
2

06

Version 3 GPU computation, PCIe communication,and MPI communication are overlapped
Figure 3 Several versions of the Cholesky
factorization algorithm

We made this decision on the basis of the preliminary evaluation, which showed that using both GPUs and CPUs slightly degraded the total performance.
We consider this to be a result of the increased use of the CPU bus and memory bandwidth, which will be investigated in detail later.

Conclusion

6This paper described a high-performance solver, SDPARA 7.5.0-G,
TABLE 1 PROBLEMS CONSIDERED IN THE EVALUATION

for large-scale, SDP problems. The key for high performance isthe acceleration of Cholesky factorization by using thousandsof GPUs. With 4,080 NVIDIA M2050 GPUs on the TSUBAME 2.0
supercomputer, our implementation achieved 533 TFlops indouble precision for a large-scale problem with m = 1,484,406.
Since SDPARA is a general solver for real problems, we not only
improved the dense matrix computation component but also
reviewed the entire software package to eliminate bottleneck
par ts that have previously been hidden. In the review, we
also made the following improvements: paralleliz ation of
the solver initialization phase (INIT) and modification of thedata structure used for the generation of Schur complement
matrix (ELEMENTS). With all these improvements, we havedemonstrated that solving SDP problems with m > 10 6 is now

Figure 4 Performance of GPU Cholesky factorization
obtained by using up to 64 nodes (192 GPUs).
Version 3 and Version 2 of CHOLESKY
are comparedpossible on modern accelerated supercomputers.
Acknowledgements
T h is r e s e a rc h w a s s u p p o r te d b y t h e J a p a n S c i e n ce a n dTechnology Agency (JST), the Core Research of Evolutionary
Science and Technology (CREST ) research projec t, and the
TSUBAME 2.0 Supercomputer Grand Challenge Program at theTokyo Institute of Technology. The authors would like to thankProf. Hayato Waki for providing us with SDP problems for QAP.
References[1]

K. Fujisawa, T. Endo, H. Sato, M. Yamashita, S. Matsuoka and
M. Nakata: High-Performance General Solver for Extremely
L a r g e -S c a l e S e m i d e f i n i t e P r o g r a m m i n g P r o b l e m s ,Figure 5 Performance of GPU Cholesky factorizationon the full TSUBAME 2.0 system (up to 1,360
nodes and 4,080 GPUs)

The proceedings of the 2012 ACM/IEEE conference on
Supercomputing, SC12 (2012)[2]

M. Yamashita, K. Fujisawa, M. Fukuda, K. Kobayashi, K.
Nakata, M. Nakata: Latest Developments in the SDPA
Family for Solving Large-Scale SDPs, M.F. Anjos and J.B.
Lasserre (eds.), Handbook on Semidefinite, Cone andPolynomial Optimization: Theory, Algorithms, Software andApplications, Springer Press, Chapter 24, 687–713 (2011)[3]

K . Fujisawa, M. Kojima, K . Nakata: E xploiting sparsit yin primal-dual interior-point methods for semidefinite
programming. Math. Prog., 79, 235–253 (1997)07

High-Performance General Solver for Extremely
Large-Scale Semidefinite Programming Problems[4]

K. Fujisawa, K. Nakata, M. Yamashita, M. Fukuda: SDPA
Project: Solving large-scale semidefinite programs. J. Oper.
Res. Soc. Japan, 50, 278–298 (2007)[5]

M. Yamashita, K. Fujisawa, M. Kojima: SDPARA: SemiDefinite
Programming Algorithm paRAllel version. Parallel Comput.,
29, 1053–1067 (2003)[6]

M. Yamashita, K. Fujisawa, M. Fukuda, K. Nakata and M.
Nakata, Parallel solver for semidefinite programming
problem having sparse Schur complement matrix, To appear
in ACM Transactions on Mathematical Software (2012)[7]

L.S. Blackford, J. Choi, A. Cleary, E. D’Azevedo, J. Demmel,
I. Dhillon, J. Dongarra, S. Hammarling, G. Henry, A. Petitet,
K . Stanley, D. Walker, R .C. Whaley: ScaL APACK Users’
Guide. Society for Industrial and Applied Mathematics,
Philadelphia (1997)[8]

P.R. Amestoy, I.S. Duff, J.-Y. L’Excellent: Multifrontal parallel
distributed symmetric and unsymmetric solvers. Comput.
Methods in Appl. Mech. Eng., 184, 501–520 (2000)[9]

I. D. Ivanov, E. de Klerk : Parallel implementation of asemidefinite programming solver based on CSDP in a
distributed memory cluster. Optim. Meth. Sof tw., 25(3),405–420 (2010)[10] S. J. Benson: Parallel computing on semidefinite programs
Preprint ANL/MCS-P939-0302 (2002)[11] M. Kojima, S. Shindoh, S. Hara: Interior-point methods
for the monotone semidefinite linear complementarity
problems. SIAM J. Optim., 7, 86–125 (1997)[12] S. Matsuoka, T. Endo, N. Maruyama, H. Sato, S. Takizawa: The
Total Picture of TSUBAME2.0. TSUBAME e-Science Journal,
GSIC, Tokyo Institute of Technology, Vol. 1, 2–4 (2010)[13] T. Endo, A. Nukada, S. Matsuoka, N. Maruyama: Linpack
Evaluation on a Supercomputer with Heterogeneous
Accelerators. Proceedings of IEEE IPDPS10, 1–8 (2010)[14] A . Petitet, R .C. Whaley, J. Dongarra, A . Clear y: HPL A
Portable Implementation of the High-Performance Linpack
Benchmark for Distributed-Memory Computers,
http://www.netlib.org/benchmark/hpl/

08

The multiple GPU calculation
for the classical spin model
Yukihiro Komura* Yutaka Okabe**Department of Physics, Tokyo Metropolitan University

With the advance of computers, large -scale and high-speed computations have become possible. In recent
years, the application of GPU to scientific computations is a hot topic in computer science. We can investigate
the phenomena such as phase transitions starting from the microscopic models in the framework of statisticalmechanics. As a numerical approach to deal with many-body problems, the Monte Carlo simulation is widely used.
Among the algorithms of Monte Carlo methods, the cluster-flip Monte Carlo method is a very efficient one, but theparallelization of the cluster-flip Monte Carlo method is not straightforward. In this paper, we describe the massive
and fast GPU computation of Swendsen-Wang multi-cluster algorithm for classical spin systems. Especially, we
implement the algorithm for multiple GPUs on the system of TSUBAME 2.0, and we realize the fast and large-scalecomputations.

Introduction

1

When we see the materials surrounding us from the microscopic
viewp oint (atomic level), they are composed of atoms or
molecules of the order of Avogadro’s number (6×1023). It is wellFigure 1 The Ising model on the two-dimensionalsquare latticeknown that water, for example, takes three states of matter,gases, liquids, and solids. The phenomenon that the statetransforms with the change of macroscopic variables, such as

W h en th e temp er at ure d ep en d en ce o f p hy si c altemperature and pressure, is called phase transition. Statisticalq u a nt i t i e s is co n s i d e r e d i n t h e f r a m e w o r k o f s t at is t i c a lmechanics is a probabilistic approach to describe the lawsmechanics, we take an average that depends on temperatureof physics among macroscopic variables, such as pressure,(Boltzmann average). For example, the total energy of thevolume, temperature, and so on, starting from the microscopic

system at a certain temperature T can be calculated throughworld. We can discuss the phase transition from the microscopic
level. When we consider the phase transition in statisticalme chanic s, it is r ather eas y to treat the magnetic phas etransition. Magnetic materials such as iron lose the properties ofmagnets when the temperature increases.
The simplest model of classical spin systems to deal
with the magnetic phase transition is the Ising model (Fig.1). On

Here, k b is the Boltzmann constant, and the summation of {S i}the lattice sites, the spin (microscopic magnet) S i is placed, andis taken over all the configurations. To take all the possible spin

the spin takes up or down direction. As a variable, S i takes ±1.

configurations means that when dealing with the Ising system

The Hamiltonian (energy) of the Ising model is given byof only 2500 spins,
= 10

752

the computational amount reaches 2 2500

; the direct calculation needs an enormous amount oftime. In the case of systems with no interaction such calculationof Boltzmann average may become simple, but systems with
Here, J (>0) is the coupling constant, and the summation is takenover the nearest neighbor pairs <i,j>. The energy is equal to -J
if the spin on the lattice site i and the nearest neighbor spin areparallel, and the energy becomes +J if the two spins are antiparallel. The important thing is that the local energy depends
on the values of the two spins, which means that there is aninteraction.

interaction cannot be simplified except for special cases.
Thus, the probabilistic Markov Chain Monte Carlosimulation has developed as a numerical technique in thefield of statistical mechanics. In the Markov Chain Monte Carloalgorithm, the method to generate states with the probabilityof the Bolt zmann distribution, Eq. (2), is called Metropolis
method [1]. If such states are generated, we can calculate thethermal average of a physical quantity at a certain temperature

09

The multiple GPU calculation
for the classical spin modelby taking a simple average. But the Metropolis method suffers

The procedure of the SW multi- cluster algorithm

from the problem of slow dynamics that it takes a long time

consists of three steps, which is shown in Fig. 2. One spin locatesto equilibrate, for example, at the temperature range nearin one mesh in Fig. 2, and the spin variables S i are expressed as S i

the phase transition point. To overcome this difficulty, several= 0, 1. There are three main steps in the SW algorithm:

algorithms, such as the cluster spin-flip algorithm [2,3] and multicanonical algorithm [4], have been proposed.

(1) Construct a bond lattice of active or non-active bonds

It goes without saying that the simulation is performedon the computer; the development of computer has changed
the role of simulation. The drastic progress of the computational
speed has made possible to treat larger system sizes, anddepending upon the temperature.
(2) The active bonds partition the spins into clusters which are
identified and labeled using a cluster-labeling algorithm.
(3) All spins in each cluster are set randomly to 0 or 1.

more realistic and more complex systems become subjects ofsimulation. At the same time, the development of numerical
techniques ranging from the computational algorithm to theanalysis method has been the matter of research all the time.
In re cent year s, the application of GPU (gr aphic
processing unit) to scientific computations has been a hot topic
in computer science. If we can efficiently use the system of
TSUMABE 2.0, a drastic advance in the field of simulation can be
expected.
As a Monte Carlo simulation of spin systems using
G PUs , Pre is e t a l . [5 ] re p o r te d t h e G PU c a l c u l at i o n o f t h eMetropolis method, and the drastic speed-up of the calculation
over the CPU calculation was realized. The Metropolis simulation
for systems with short-range interaction is a problem of local
computation, and the parallelization is rather simple. Instead,the cluster-flip method is a non-local computation, and theparallelization is not straightforward. We have succeeded inrealizing the GPU-based calculation for cluster-flip method;
fast and large-scale computations of the cluster-flip Monte
Carlo method are highly desired. In this paper, we explain the

Figure 2 The procedure of the SW multi-cluster algorithm

calculation of cluster-flip Monte Carlo method using a single
GPU and that using multiple GPUs, and the performance of ouralgorithm on the TSUBAME 2.0 is reported.

Cluster-flip Monte Carlo algorithm

is regarded as one Monte Carlo Step. The Hoshen-Kopelman

2

Let us s t ar t w ith e xp laining the clus ter- f lip M onte C ar lo
method. This algorithm is an algorithm that the clusters of
the spins are generated and then the spins in each cluster
are updated at once. The cluster-flip algorithm successfully
reduced the problem of slow dynamics. In this paper, we deal
with the multi-cluster-flip algorithm proposed by Swendsen
and Wang (SW) [2] ; they applied the Fortuin and Kasteleyn [6]representation to identify clusters of spins.

10

The procedure from the step (1) to the step (3) is repeated, whichalgorithm [7], which is a special version of the class of union-andfind algorithms [8], is often used for an efficient cluster-labelingalgorithm. Integer labels are assigned to each spin in a cluster.
Each cluster has its own distinct set of labels. The proper labelof a cluster, which is defined to be the smallest label of any spinin the cluster, is found by the following function. The array label
is used, and if label is a label belonging to a cluster, label[label]is the index of another label in the same cluster which has asmaller value if such a smaller value exists. The proper label forthe cluster is found by evaluating label[label] repeatedly until
the label[label] is unchanged.

The calculation with a single GPU

3

Now, we compare the p er formance on the t wo dimensional Ising model with CPU to that with a single GPU. We
have tested the performance of our code on NVIDIA GeForce
GTX580. For comparison, we run the code on a current CPU,Intel(R) Xeon(R) CPU W3680 @ 3.33 GHz; we use only one core of

When we extend the CPU calculation to the GPU calculation, we

the CPU. For compiler, we have used gcc 4.1.2 with option -O3.

should keep in mind to change the algorithm from a sequential

We employ the HK algorithm as the method of cluster labelingcomputation to a parallel computation. This shift is very importantwith CPU. Since the cluster size depends on the temperature,

to realize the GPU calculation. In the SW multi-cluster algorithm

we compare the performance at the critical temperature of thewith GPU, the procedure from the step (1) to the step (3) remains

two-dimensional Ising model k bTc /J= 2.269... [12]. We show theunchanged. Thus, we check whether the SW multi- clustercomputational time per one Monte Carlo Step (millisecond) inalgorithm with CPU is directly applicable to that with GPU. In the

Table 1. There, the linear system sizes L are L = 256, 1024 andstep (1), each bond can be generated independently. In the step4096. We can see from Table 1 that the computational times with(2), the method of cluster labeling with CPU cannot directly be

GPU are highly accelerated compared to that with CPU. We also

applied to the single GPU calculation since the HK algorithm is aconfirm that the computational speed using the algorithm ofsequential computation. In the step (3), each spin can be updated

Kalentev et al. is superior to that of Hawick et al. for all system sizes.

independently if the new spin states are given in advance. From
the above, the procedure with CPU in the steps (1) and (3) aredirectly applicable to that with GPU. However, we should modify
the algorithm of cluster labeling in the step (2) so as not to use

The calculation with multiple GPUs

the sequential computation.

4

The cluster labeling algorithm with a single GPU
was proposed by Hawick et. al.[9] and the refinement of that by

The system size we can treat is limited for a single GPU. If we use

Kalentev et al. [10] Those cluster algorithms are named the label

a multiple GPU system such as TSBAME 2.0, a very large system

equivalence algorithm. We employ the two label equivalence

can be treated directly, and the fast computation becomesalgorithms for the method of cluster labeling with a single

possible. The GPU calculation itself is a parallel calculation, but

GPU [11] . The concept of the label equivalence algorithm isthere are several problems that we should bear in mind whenbasically equal to the HK algorithm; that is, each thread updates

performing multiple GPU calculation. First, the communication

the label[label] until the label[label] is unchanged. However,between multiple GPUs should be considered. Second, thethe label equivalence algorithm repeats the kernel call withcalculation with a single GPU is the shared memory calculation,

the label updates, which is different from the case of the HK

whereas the calculation with multiple GPUs is the distributed

algorithm. We note that the label equivalence algorithms arememory calculation. The difference of memory type becomes a

explained in Ref. [11] in detail.

serious problem in the extension of SW multi-cluster algorithm
from a single GPU to multiple GPUs.
As a model, we take the two-dimensional Ising model
again. The data transfers on three levels are needed, that is, thedata transfer from GPU to CPU, that between CPUs and thatfrom CPU to GPU. We use API of CUDA for the communication
between CPU and GPU, and MPI library for the communication
between CPUs. We arrange the total lattice with a super-lattice
structure using multiple GPUs. Each GPU has the informationof spins on a sub-lattice together with the arrays to preserve
the data of surrounding boundary layers and to transfer the

Table 1 The computational times per one Monte Carlo
Step (millisecond) for the two-dimensional Isingmodel at the critical temperature with a CPU anda single GPU.

data of boundary layers. We illustrate the case of a 2 × 2 superlattice structure using 4 GPUs in Fig. 3. To arrange the data ofsurrounding boundary layers in each GPU, the calculation in
the step (1) can be executed without the data transfer with

11

The multiple GPU calculation
for the classical spin model

its neighbors. In the step (2), we employ a two-stage process

The cluster labeling algorithm with multiple CPUs isof cluster labeling. After the cluster labeling within each GPU

realized by using one label variable. On the other hand, we use

is finished, we check the bond between the sites of ‘‘nearest-

two label variables, that is, the local label and the global label,neighbor ’’ GPU. In the inter- GPU labeling step, we should

to realize the cluster labeling algorithm with multiple GPUs. To

consider the distributed memory, which will be discussed next.

resolve the problem in the step (3), we store the new spin state

The treatment for distributed memory should be consideredin the label variable by using the bit manipulation.

also in the step (3).

We have tested the performance of the multiple GPU
calculation of the SW cluster-flip algorithm on TSUBAME 2.0 withNVIDIA Tesla M2050 GPUs using CUDA 4.0 and openMPI 1.4.2.
As a model, we use the two-dimensional Ising model at thecritical temperature. We give the double logarithmic plot of theaverage amount of spin flips per nano second as a function of
the number of GPUs in Fig. 4. For the size of sub-lattice of eachGPU, not only the data for 4096 × 4096 but also those for 2048
× 2048 and 1024 × 1024 are plotted. The dependences on thenumber of GPUs with fixing the total linear system size as L =4096, 8192, 16384, 32768 and 65536 are shown. In the inset, thedependence on the number of GPUs with fixing the sub-lattice

Figure 3 The 2 × 2 super-lattice structure for 4 GPUs
and the information on a single GPU.

size as 4096 × 4096 is shown. We see that the performanceof our code increases as a power of the number of GPUs. The
coefficient of the power is estimated as 0.91 by the best-fitted
curve shown in the inset. Since it takes extra time to diffuse

If we directly apply the SW multi-cluster algorithm

the label number from one GPU to all GPUs, the perfect weakwith a single GPU to that with multiple GPUs, the problems

scalability, that is, the power of 1.0, is not achieved, but we canoccur in the steps (2) and (3). In the step (2) of cluster labelingsay that the weak scalability is well satisfied for our code.

with a single GPU, each thread updates the label[label] until
the label[label] is unchanged. On the other hand, in the step(2) with multiple GPUs, each thread has a possibility to accessthe array on different GPU memory. And, this access to different
GPU memory generates errors in the step (2) because each GPU
cannot get access to the different GPU memory. If we directly
apply the step (3) of spin flip with a single GPU to that withmultiple GPUs, each GPU should check the label number in eachGPU after the step (2) and make a decision about the targets
of communication in all GPUs; then the new spin states aretransferred from those targets to each GPU. But, that calculation
with multiple GPUs needs an enormous amount of costs.
Some at tempts for the cluster labeling algorithmwith multiple CPUs have been proposed: the master-slaves
method [13], the method to prepare a global label array [14], andthe method to exchange the boundary information with its
neighbors until the all labels are unchanged (named a relaxation
method) [15]. In a relaxation method, there are some versions; themethod to set up the local root [15] and the method to generate
the area of overlap [16]. In this paper, we employ the relaxation
method since the amount of memory on a single GPU is not so
much.

12

Figure 4 Plot of spin flips per a nano second withchanging the number of GPUs for the twodimensional Ising model at the criticaltemperature on the TSUBAME 2.0.

References

5

Conclusion[1] N. Metropolis, A.W. Rosenbluth, M.N. Rosenbluth and A.H. Teller, E. Teller:
Equation of State Calculations by Fast Computing Machines, J. Chem. Phys.
Vol. 21, pp. 1087-1092 (1953)

We proposed the GPU-based calculation of the SW multicluster algorithm, which is one of the Markov Chain Monte Carloalgorithms in statistical mechanics, with a single GPU [11] andmultiple GPUs [17]. The cluster labeling algorithm with a CPU is not[2] R.H. Swendsen and J.S. Wang: Nonuniversal critical dynamics in Monte Carlosimulations, Phys. Rev. Lett. Vol. 58, pp. 86-88 (1987)[3] U. Wolff: Collective Monte Carlo Updating for Spin Systems, Phys. Rev. Lett.
Vol. 62, pp.361-364 (1989)

directly applicable to that with GPU. However, we have realized[4] F. Wang and D.P. Landau: Phys. Rev. Lett. Vol. 86 pp. 2050-2053 (2001).

the SW cluster algorithm by using the two label equivalence[5] T. Preis, P. Virnau, W. Paul and J.J. Schneider: GPU accelerated Monte Carloalgorithms proposed by Hawick et al.

[9]and that by Kalentev et al.

[10]

.

Moreover, we have extended the SW cluster algorithmwith a single GPU to that with multiple GPUs. We cannot directly

simulation of the 2D and 3D Ising model, J. Comp. Phys. 228 pp. 44684477(2009).
[6] C. M. For tuin and P. W. Kasteleyn: On the random-cluster model: I.

apply the cluster labeling algorithm with a single GPU to that with

Introduction and relation to other models, Physica Vol 57, pp. 536-564 (1972)

multiple GPUs because the memory types of a single GPU and[7] J. Hoshen, and R. Kopelman: Percolation and cluster distribution. I. Cluster

multiple GPUs are different. By using the relaxation method andmultiple labeling technique and critical concentration algorithm, Phys. Rev.

two label variables, we have realized the cluster labeling algorithmwith multiple GPUs in the SW multi-cluster algorithm. In the stepof cluster labeling with multiple GPUs, we employ a two-stage

B. Vol. 14, pp.3438-3445 (1976).
[8] T. H. Cormen, C. E. Leiserson, R. L. Rivest and C. Stein: Introduction to
Algorithms, MIT Press, third edition, 2009.

process; that is, we first make the cluster labeling within each[9] K.A. Hawick, A. Leist and D. P. Playne: Parallel Graph Component Labelling

GPU, and then the inter-GPU labeling is performed with referencewith GPUs and CUDA, Parallel Computing Vol. 36, pp. 655-678 (2010)

to the label on the first-stage. We use some tricks for the step of[10] O. Kalentev, A. Rai, S. Kemnitzb and R. Schneider: Connected component

spin update; we include the information on new spin state in thelabeling on a 2D grid using CUDA, J. Parallel Distrib. Comput. Vol. 71, pp.

process of cluster labeling. We have tested the performance of the615-620 (2011)

SW multi-cluster algorithm with multiple GPUs on TSUBAME 2.0 for[11] Y. Komura and Y. Okabe: GPU-based Swendsen-Wang multi-cluster

the two-dimensional Ising model at the critical temperature. We

algorithm for the simulation of two-dimensional classical spin systems,

have shown a performance with good scalability for multiple GPU
computations although the cluster-flip algorithm is a nonlocal
problem which is not so easy for parallelization.

Comp. Phys. Comm. Vol. 183, pp. 1155-1161 (2012)[12] L. Onsager: Crystal Statistics. I. A Two-Dimensional Model with an OrderDisorder Transition, Phys. Rev. Vol. 65, pp. 117-149 (1944)

The SW cluster-flip algorithm using multiple GPUs can[13] M. Bauerfeind, R. Hackl, H.G. Matuttis, J. Singer, I. Morgenstern, 3D Ising

be applied to a wide range of spin problems. As an example,

model with Swendsen–Wang dynamics: a parallel approach, Physica A Vol.

we have already studied the large-scale calculation for the212, pp. 277–298 (1994)

two-dimensional classical XY model[18]. Moreover, the cluster[14] J.M. Constantin, M.W. Berry, B.T. Vander Zanden, Parallelization of thelabeling using multiple GPUs can be applied to several problems

Hoshen–Kopelman algorithm using a finite state machine,Internat. J.

including the percolation, the image processing, and so on.

Supercomput. Appl. High Perform. Comput. Vol. 11, pp. 34–48 (1997)

We are now attempting to refine the cluster labeling algorithm[15] M. Flanigan, P. Tamayo, Parallel cluster labeling for large-scale Monte Carlowith multiple GPUs and to develop the implementation of thealgorithms on the problems in various fields.

simulations, Physica A Vol. 215, pp. 461–480 (1995)[16] J.M. Teuler and J.C. Gimel :A direct parallel implementation of the Hoshen–
Kopelman algorithm for distributed memory rchitectures Comp. Phys.

Acknowledgements

Comm. Vol. 130, pp. 118-129 (2000)[17] Y. Komura and Y. Okabe: Multi-GPU-based Swendsen–Wang multi-cluster

We thank Takayuki Aoki, Yusuke Tomita, and Wolfhard Janke foralgorithm for the simulation of two-dimensional q-state Potts model, Comp.

valuable discussions. This work was supported by a Grant-in-Aid

Phys. Comm. Vol. 184, pp. 40-44 (2012)

for Scientific Research from the Japan Society for the Promotion[18] Y. Komura and Y. Okabe: Large-Scale Monte Carlo Simulation of Two-of Science. The computations have been per formed on the

Dimensional Classical XY Model Using Multiple GPUs, J. Phys. Soc. Jpn., Vol.81,

TSUBAME 2.0 system of the Tokyo Institute of Technology.

113001 (2012).-13

Knowledge Acquisition froma Large Web Corpus andits Applications
Daisuke Kawahara* Sadao Kurohashi**Graduate School of Informatics, Kyoto University

In Science fiction novels and movies, robots communicate naturally with humans. In order to realize such robots,
it is necessary to make them understand language and also to provide them with knowledge that is common tohumans. However, it is very difficult to express such knowledge manually, owing to the vastness thereof. As such,the development of robots has been affected by this deadlock for a long time. Owing to the recent progress incomputer networks, and the web in particular, a large number of texts has been accumulated and it becomes
possible to acquire knowledge from such texts. This article introduces automatic knowledge acquisition from a large
web text collection and linguistic analysis based on the acquired knowledge.

Introduction

1

As the acquired linguistic knowledge, this ar ticle
describes case frames, which represent aggregated predicateargument structures, namely "who does what." This is one of basic
linguistic knowledge known by humans. We compile case frames
in the following way[1]:R e centl y, we have b e co m e accus to m e d to usin g sp e e chinterfaces for smart phones, whitch facilitate the execution of

1. Extraction of Japanese sentences from a web pagecertain specific tasks in practice. However, it is still not possible

collection (corpus creation)

for computers to communicate naturally with humans. To realize

2. Syntactic parsing of the corpus

such computers, it is necessary to develop computer systems that

3. Filtering and clustering of syntactic parses to produce

can understand and use natural language at the same level ascase frames

humans.
The following two sentences are classical examples ofsyntactic ambiguities in Japanese sentences.

This time, we transmitted the corpus that had been created onother clusters to TSUBAME, and performed the processes ofsyntactic analysis of the corpus and case frame compilation on（1）a . クロールで泳ぐ少女を見た
b. 望遠鏡で泳ぐ少女を見た

TSUBAME.
In the remainder of this article, we first describe themethod of case frame compilation. Then, we explain linguistic

These two sentences have the same style, but "クロールで"analysis based on acquired case frames. Finally, we report the(crawl) and "望遠鏡で" (telescope) have different modifying

experimental results on TSUBAME.

heads (the underlined part modies the doubly-underlined part).
This kind of syntactic ambiguity cannot be resolved by Japanesegrammar, but by precisely considering the relations between

Automatic Compilation ofwords. To do that, linguistic knowledge of words such as "クローCase Frames

ルで泳ぐ" (swim crawl) and "望遠鏡で見る" (see with telescope) is2required. However, the amount of knowledge required is too vast
to be expressed manually.
In recent years, a massive number of tex ts can beobtained from the web, and studies that automatically acquire

A case frame represents the relations between a predicate andits arguments for each usage of the predicate. For example, thefollowing is a case frame for the predicate " 焼く " :linguistic knowledge from such a large text collection havearisen. We automatically acquire the linguistic knowledge froma huge tex t collection crawled from the web, and develop（2）{ 私、人、… ｝が{ パン、肉、… ｝を{ オーブン、フライパン、… ｝で焼く

an analyzer based on the acquired knowledge. This ar ticle

14

introduces that this knowledge acquisition process was carried

S u c h c a s e f r a m e s a r e a u to m a t i c a l l y a cq u i r e d b y

out in an extremely short period using the large computational

applying s yntac tic parsing to a huge corpus, ex trac tingresources of TSUBAME.

syntactically unambiguous relations between a predicate andits arguments from the resulting parses and clustering them.
Syntactically unambiguous relations consist of those arguments

Predicate-argument Structure Analysis

that have only one modif ying candidate head based on thebased on Case Frames

gr ammatical rules of the s y ntac tic p ar ser. In the following

3

examples, the underlined arguments demarcated with a ○ areextracted as unambiguous ones:

By integrating the automatically compiled case frames into
syntactic parsing, we can now resolve the syntactic ambiguities（3）a. 今日は石窯で ○ パンを ○ 焼いています。

in examples (1a) and (1b)[2]. For example, it is difficult to acquireb. 80 種類ものパンを ○ 焼いていますが、…

the relation " クロ ールで 泳ぐ " (swim crawl), which is required toc. その後、パンを × 焼いた余熱を利用し…

correctly analyze example (1a), from a newswire corpus comprisingd. 直径が × 15センチのケーキを焼きます

tens of thousands of sentences. However, it is included in thecase frames compiled from a large corpus. That is to say, wide-

In ( 3 a) and ( 3 b), since the sentence end and the subordinate

coverage knowledge on relations between content words is

clause " ～ が " are considered to be strong boundaries, thesuccessfully captured in the case frames. These case frames enablearguments " 石 窯でパンを" (masonry oven, bread) and "パンを"to disambiguate syntactic ambiguities as in example (1).

(bread), which unambiguously modify " 焼く" (bake), are extracted
from these sentences, respectively.

Further, it is possible to detect the hidden cases of topicmarked phrases and clausal modifiees. For example, in the two

On the other hand, the underlined arguments

sentences in (1), "少女 " (girl) has a "が" (nominative) relation to " 泳demarcated with a × in (3c) and (3d) are not extracted due to theぐ" (swim). This relation can be analyzed by predicate-argument

ambiguity of modifying heads. In (3c), the head of the underlined

structure analysis based on the case frames.

argument "パンを" (bread) is correctly analyzed as "焼いた" (baked),

This kind of method, which learns relations between

but this argument is not extracted because it has two candidate

words from a large corpus and exploits them in analysis, has beenheads " 焼いた" (baked) and " 利用し" (use). In (3d), the underlined

applied to other languages [4, 7].

argument " 直径が" (diameter) is not extracted either, because it
has two candidate heads " 15センチ" (15cm) and " 焼きます" (bake).
Since the syntactic parser incorrectly analyzes the head of " 直径が" (diameter) as " 焼きます" (bake), extracting such ambiguous
relations would produce erroneous case frames. This filtering canbe done on the basis of the grammatical rules in the syntactic
parser, which handle such exceptional phrases.
The biggest problem in automatic case framecompilation is sense ambiguity of predicates. In case of " 焼く," itis necessary to distinguish expressions with different meanings
and usages such as "パン を 焼 く" (bake bread) and " 手 を 焼 く"
(have difficulty). To deal with this problem, predicate-argument
structures are distinguished by coupling a verb and its closest
argument. This is because the closest argument is thought tostrongly influence the meaning and usage of the predicate. For" 焼 く," by making a couple of "パン を 焼 く" (bake bread), " 肉 を焼く" (broil meat) and " 手 を 焼く" (have difficulty), it is possible

Table 1

Acquired case frames for yaku. Example
words are expressed only in English due
to space limitation. The number following
each word denotes its frequency. `CS'
indicates a case slot.

to separate predicate-argument structures according to their
meanings and usages and to acquire a case frame for each ofthem. Further, by applying clustering based on word similarities,
similar case frames such as "パンを焼く" (bake bread) , and " 肉を焼く" (broil meat) can be merged to produce final case frames asgiven in Table 1. The word similarities are measured on the basisof distributional similarities calculated from the huge corpus[3].

15Knowledge Acquisition froma Large Web Corpus and its Applications

Experimental Results

4First, we extracted Japanese sentences from three billion web
pages, eliminated duplicated sentences and obtained a corpus
comprising approximately 15 billion sentences. This extraction
had been conduc ted on other PC clusters and the ex trac ted
corpus was transmitted to TSUBAME. We then applied syntactic
parsing to this corpus using TSUBAME. To carry out this process,we needed computational power of approximately 200,000 CPU
core hours. The smallest unit in this process is a sentence, andthus we were able to implement it in an embarrassingly parallel
way by split ting the corpus into small processing units eachcousisting of approximately 10,000 sentences.
Secondly, we acquired case frames from the parses of
15 billion sentences using the method of case frame compilation
described in section 2.1 As a result, we acquired case frames forapproximately 40,0 0 0 predicates. To perform this process, we
needed computational power of approximately 10 0, 0 0 0 CPU
core hours. By integrating the acquired case frames into our base

Figure 1 Changes in accuracy and coveragefor varying sizes of corpora

syntactic parser, we developed a predicate-argument structure
analyzer.
Thirdly, to investigate the relationship between corpus
size and the accuracy of the analyzers, we created a corpus of 1.5
million, 6 million, 25 million, 100 million, 400 million, 1.6 billion
and 6.4 billion sentences by sampling the above corpus consistingof 15 billion sentences. We conducted case frame compilationfrom each corpus and developed analyzers based on these caseframes. Figure 1 shows the changes in accuracy and coveragefor the varying sizes of corpora. In this figure, "coverage of caseframes" represents how many predicate-argument structuresin the test sentences are included in the case frames, "synonymrecognition" means the recognition of synonymous phrases suchas "景気が悪化する" (economy gets worse) and "景気が冷え込む"(economy gets cold), and "anaphora resolution" means an analysis
of detec ting omit ted nouns such as the subjec t of a verb[ 5 ].
From this figure, we can see the improvements in accuracy andcoverage as the size of corpus increases.
To conduct the above computation on TSUBAME, weused the grid shell, GXP2 , which has been developed at the Tauralaboratory at the University of Tokyo. As the queue for TSUBAME,
we mainly used the reservation queue (H queue) with at most 300
nodes (3,600 CPU cores) in parallel.

16

1

The acquired case frames can be searched from http://nlp.ist.i.kyoto-u.ac.jp/index.php? 京 都 大 学 格フレーム.

2

http://www.logos.ic.i.u-tokyo.ac.jp/gxp/

Conclusion

5

References[1]

Daisuke K awahara and Sadao Kurohashi. Case f rame
co m p i l a t i o n f r o m t h e w e b u s i n g h i g h p e r f o r m a n ce
computing. In Proceedings of LREC2006, 2006.

In this article, we have described case frame compilation from

[2]

Daisuke Kawahara and Sadao Kurohashi. A fullylexicalized

a huge tex t collec tion collec ted from the web and syntac tic

probabilis tic mo del for Japanese s yntac tic and casea n d p r e d i c a te - a r g u m e nt s t r u c t u r e a n a l y s is b a s e d o n t h e

structure analysis. In Proceedings of HLT-NAACL2006, pp.

acquired case frames. The process of case frame compilation
was per formed in an ex tremely shor t period using the large

176-183, 2006.
[3]

computational resources of TSUBAME.
A s a n a p p l i c a t i o n o f t h e a u to m a t i c a l l y a cq u i r e dDekang Lin. Automatic retrieval and clustering of similarwords. In Proceedings of COLINGACL98, pp. 768-774, 1998.

[4]

David McClosk y, Eugene Charniak, and Mark Johnson.

linguistic knowledge and the analyzer based on this knowledge,

Effective self-training for parsing. In Proceedings of HLT-

we have been developing the search engine infrastruc ture,

NAACL2006, pp. 152-159, 2006.

TSUBAKI[6]. The biggest characteristic of TSUBAKI is the precise[5]

Ryohei Sasano and Sadao Kurohashi. A discriminative

search exploiting synonym knowledge and linguistic structure

approach to japanese zero anaphora resolution with large-(as in Figure 2). In this way, the automatically acquired linguistic

scale lexicalized case frames. In Proceedings of IJCNLP2011,

k nowledge and deep analysis can be exploited in prac ticalapplications effectively.

pp. 758-766, 2011.
[6]

Keiji Shinz ato, Tomohide Shibata, Daisuke Kawahara,

The accuracy of the analyses and the coverage of caseand Sadao Kurohashi. Tsubaki: An open search engine

frames are expected to be improved with an increasing corpus

infrastruc ture for developing information accesssize. In future, we will tr y to increase the size of the corpus tomethodology. Journal of Information Processing, Vol. 20, No.

acquire more knowledge. For this purpose, we will use TSUBAME
on a larger scale.

1, pp. 216-227, 2012.
[7]

Jun Suzuki, Hideki Isozaki, Xavier Carreras, and Michael
Collins. An empirical study of semi-supervised structured
conditional models for dependency parsing. In Proceedingsof EMNLP2009, pp. 551-560, 2009.

● TSUBAME e-Science Journal No.7
Figure 2 A search result of search engine TSUBAKI

Published 26/12/2012 by GSIC, Tokyo Institute of Technology ©
ISSN 2185-6028
Design & Layout: Kick and Punch
Editor:TSUBAME e-Science Journal - Editorial room
Takayuki AOKI, Thirapong PIPATPONGSA,
Toshio WATANABE, Atsushi SASAKI, Eri Nakagawa
Address: 2-12-1-E2-1 O-okayama, Meguro-ku, Tokyo 152-8550
Tel:
E-mail:URL:+81-3-5734-2087 Fax: +81-3-5734-3198
tsubame_j@sim.gsic.titech.ac.jp
http://www.gsic.titech.ac.jp/

17International Research Collaboration
The high performance of supercomputer TSUBAME has been extended to theinternational arena. We promote international research collaborations using
TSUBAME between researchers of Tokyo Institute of Technology and overseas
research institutions as well as research groups worldwide.

Recent research collaborations using TSUBAME
1. Simulation of Tsunamis Generated by Earthquakes using Parallel
Computing Technique
2. Numerical Simulation of Energy Conversion with MHD Plasma-fluid Flow
3. GPU computing for Computational Fluid Dynamics

Application Guidance
Candidates to initiate research collaborations are expected to conclude
MOU (Memorandum of Understanding) with the partner organizations/departments. Committee reviews the “Agreement for Collaboration” for joint
research to ensure that the proposed research meet academic qualifications
and contributions to international society. Overseas users must observe
rules and regulations on using TSUBAME. User fees are paid by Tokyo Tech’s
researcher as part of research collaboration. The results of joint research are
expected to be released for academic publication.

Inquiry
Please see the following website for more details.
http://www.gsic.titech.ac.jp/en/InternationalCollaboration