量子アニーリングと機械学習の交差点

2021/11/12

門脇正史(株)デンソー AI研究部 データサイエンス研究室
量子コンピューティング研究課

量子アニーリング (Quantum Annealing) ?
𝐽𝑖𝑗 ∈ ℝ, 𝑥𝑖 ∈ {−1,1}

二値変数のコスト関数ℋ𝐶 = − σ𝑖<𝑗 𝐽𝑖𝑗 𝑥𝑖 𝑥𝑗𝜎𝑖𝑧 , 𝜎𝑖𝑥 : パウリ行列

を、量子力学に拡張 (𝑥𝑖 → 𝜎𝑖𝑧 )し、横磁場(~量子ゆらぎ)項 σ𝑖 𝜎𝑖𝑥 を加える。
(横磁場からコスト関数への変形をパラメータ化したい)
1𝑧𝑧𝑥ℋ𝑄 𝑡 = −𝐴 𝑡 σ𝑖<𝑗 𝐽𝑖𝑗 𝜎𝑖 𝜎𝑗 − 𝐵 𝑡 σ𝑖 𝜎𝑖𝐵 𝑡

𝐴 𝑡

初期状態ȁ𝜓 0 ۧを

ℋ𝑄 0 = σ𝑖 𝜎𝑖𝑥0𝑡の基底状態*にえらび、シュレディンガー方程式によって系を時間発展する。
𝑑𝑖 𝑑𝑡 ȁ𝜓 𝑡 ۧ = ℋ ȁ𝜓 𝑡 ۧ

断熱的に(~とってもゆっくりと)時間発展すると、断熱定理によりコスト関数の最適解が求まる。
つまり、自明なハミルトニアンの解(すべての状態の重ね合わせ)から、非自明なハミルトニアンの解(知りたいコスト関数の最適解)が求まる。
ちなみにシミュレーテッドアニーリングは、自明な分布(均一分布)から、非自明な分布(最適解でのデルタ関数)が求まる。

*基底状態(最小固有値の固有状態)はすべての計算基底の重ね合わせになっている(~全解空間の探索)

量子アニーリングと機械学習の交差点 / 2021/11/12 @ IBIS2021© DENSO CORPORATION All Rights Reserved.

2 / 39量子アニーリングで何を目指す？
製品

データサイエンス

機械学習(予測)

最適化(制御)AI
モノづくり
設計/開発/生産
・素材・部品・生産・物流

サービス

モビリティ・配車・配送【量子アニーリング業界】
1. まだまだ基礎研究• もっと最適化を頑張る: 新しい相互作用、新しい制御、ダイナミクス• 最適化以外の出口: 量子多体系シミュレーション
2. 離散最適化ソルバーとして使う (イジングソルバーも含めて)
量子アニーリングと機械学習の交差点 / 2021/11/12 @ IBIS2021© DENSO CORPORATION All Rights Reserved.

今日の話

3 / 39Future of product design
Data generation••••

Model generation

experiments
simulations
robotics
database

VQAs

Variational
Quantum
Algorithm

Robot Scientist
Cyber Physical System• machine learning

QML
Quantum-powered
design automation

Hypothesis generation• optimization of cost functions
量子アニーリングと機械学習の交差点 / 2021/11/12 @ IBIS2021© DENSO CORPORATION All Rights Reserved.

QuantumMachine
Learning

QA

QuantumAnnealing
4 / 39

1. ブラックボックス最適化••

非可逆行列圧縮
基板設計最適化

2. 量子アニーリング•

機械学習 > 量子 な話

量子 > 機械学習 な話

次世代量子アニーラ向けアルゴリズム

量子アニーリングと機械学習の交差点 / 2021/11/12 @ IBIS2021© DENSO CORPORATION All Rights Reserved.

5 / 39Black-box optimization (BBO)• Black-box optimization (BBO) of is a method to optimize unknown costfunctions.

• Bayesian optimization (BO) is a well-known BBO algorithm, that
acquire data points sequentially based on a model constructed frompre-acquired data set.
• This is a Bayesian approach of experimental design.

• BO iterates the 3 processes;•••

input-output data from difficult-to-evaluate function, such as experiments
construct a regression model (e.g. Gaussian process) & an acquisition
function by using the regression model
find a new data point that maximizes the acquisition function• BO is designed for “continuous” explanatory variables.
• BBO for “integer” variables (combinatorial optimization) is not studied
so far.

量子アニーリングと機械学習の交差点 / 2021/11/12 @ IBIS2021© DENSO CORPORATION All Rights Reserved.

6 / 39FMQA algorithm

Kitai et al (Phys Rev Res 2020)
10.1103/PhysRevResearch.2.013319

Maximizeacquisition
function

Experimental(Input-output) dataModelingacquisition
function

Parameter for sparseness
量子アニーリングと機械学習の交差点 / 2021/11/12 @ IBIS2021© DENSO CORPORATION All Rights Reserved.

7 / 39

Nanostructure optimization of metamaterial

Kitai et al (Phys Rev Res 2020)
10.1103/PhysRevResearch.2.013319

The goal is to design complex structures of wavelength selective radiatorswith the thermal atmospheric transparency window.

量子アニーリングと機械学習の交差点 / 2021/11/12 @ IBIS2021© DENSO CORPORATION All Rights Reserved.

8 / 39Bayesian Optimization of CombinatorialStructures (BOCS)
Baptista and Poloczek (ICML2018)arXiv:1806.08838

5. Acquisition function
6. Maximization
7. Input-output data

sparsity-inducing prior(horseshoe prior)
量子アニーリングと機械学習の交差点 / 2021/11/12 @ IBIS2021© DENSO CORPORATION All Rights Reserved.

9 / 39

BOCS on random Binary Quadratic Programming
instances
Baptista and Poloczek (ICML2018)arXiv:1806.08838

Other algorithms
EI: Expected Improvement
PS: Sequential Monte Carlo Particle Search
MLE: Maximum Likelihood Estimate
OLS: Oblivious Local Search
SA: Simulated Annealing
SMAC: Sequential Model-based Algorithm Configuration

BOCS algorithms
BOCS-SDPBOCS-SA

量子アニーリングと機械学習の交差点 / 2021/11/12 @ IBIS2021© DENSO CORPORATION All Rights Reserved.

10 / 39

BOCS-QA

Koshikawa et al (J Phys Soc Jpn 2021)
10.7566/JPSJ.90.064001

Finding ground state of diluted Sherrington-Kirkpatrick model fromsequentially obtained state-energy data set without knowing the Hamiltonian.
1ℋ = − ෍ 𝐽𝑖𝑗 𝜎𝑖 𝜎𝑗
𝐽𝑖𝑗 ~ 1 − 𝜌 𝒩 0, +0 + 𝜌𝒩 0,1𝑁𝑖<𝑗better solution

N=20, 50 instances x 10 runs

量子アニーリングと機械学習の交差点 / 2021/11/12 @ IBIS2021© DENSO CORPORATION All Rights Reserved.

11 / 39

FMQA x photonic crystal

量子アニーリングと機械学習の交差点 / 2021/11/12 @ IBIS2021© DENSO CORPORATION All Rights Reserved.

12 / 39

Integer decomposition (ID) 整数基底分解法 (非可逆圧縮)
Decompose a matrix into an integer matrix and a real matrix

𝐷𝑁𝑊

𝐷𝐾

≈

𝑉(𝑀, 𝐶)

Ambai and Sato (ECCV2014)= 𝑁 𝑀

𝐶argmin 𝑊 − 𝑉 2𝑀∈ −1,1 𝑁𝐾
𝐶∈ℝ𝐾𝐷

𝑊: target weight matrix
𝑉: approximated matrix (function of 𝑀 and 𝐶)
𝑀: integer matrix 𝑀 ∈ −1,1 𝑁𝐾
𝐶: real matrix 𝑐 ∈ ℝ𝐾𝐷
The real matrix 𝐶 is a dominant part of the memory footprint.
Memory size reduction 𝑁 × 𝐷 → 𝐾 × 𝐷 & speed up in matrix operations.
However, finding 𝑀 and 𝐶 is a mixed integer non-linear programming(MINLP) and NP hard problem.
量子アニーリングと機械学習の交差点 / 2021/11/12 @ IBIS2021© DENSO CORPORATION All Rights Reserved.

13 / 39

Implementation of ID

Ambai and Sato (ECCV2014)

Repeated 1-rank approximation(=minimizing residuals)
𝒄1𝐶𝑀

= 𝒎1

𝒄2+ 𝒎2

1-rank approx.
of 𝑀𝐶

量子アニーリングと機械学習の交差点 / 2021/11/12 @ IBIS2021© DENSO CORPORATION All Rights Reserved.

𝐾𝒄3+ 𝒎3+ …

1-rank approx.
of 𝑀𝐶 − 𝒎1 𝒄1

1-rank approx.
of 𝑀𝐶 − 𝒎1 𝒄1 − 𝒎2 𝒄2

Residual

Residual

14 / 39

Applications of ID
1. Image recognition•Linear SVM kernels
Mitsuru Ambai and Ikuro Sato, SPADE: Scalar Product Accelerator byInteger Decomposition for Object Detection (ECCV2014)

2. Voice recognition•Deep neural networks
太刀岡勇気, 安倍満, 整数基底分解法と量子化による音響モデルの圧縮 (日本音響学会2018年春季研究発表会)Sound
input・・・

Prediction
output= 𝑓 𝑊𝑇𝑥𝑓 : activation function

ID demonstrated 1/20 memory footprint,
x15 acceleration with 1.5% increase of the error.

量子アニーリングと機械学習の交差点 / 2021/11/12 @ IBIS2021© DENSO CORPORATION All Rights Reserved.

15 / 39

Collaboration with M. Ambai

BBO for ID
We can convert the MINLP toa nonlinear integer programming (NLIP).
In ID, for given 𝑀,

𝐶 = 𝑀𝑇 𝑀 −1 𝑀𝑇 𝑊
𝑉 𝑀, 𝐶 = 𝑀𝐶 = 𝑀 𝑀𝑇 𝑀 −1 𝑀𝑇 𝑊 = 𝑉 𝑀
The problem is written asargmin 𝑊 − 𝑉 2

𝑀∈ −1,1 𝑁𝐾= argmin 𝑊 − 𝑀 𝑀𝑇 𝑀 −1 𝑀𝑇 𝑊 2
𝑀∈ −1,1 𝑁𝐾

MINLPの一般的な解法としてはBenders分解法などがある。

Non-linear & non-QUBO combinatorial optimization problem量子アニーリングと機械学習の交差点 / 2021/11/12 @ IBIS2021© DENSO CORPORATION All Rights Reserved.

16 / 39

How to solve it?

argmin 𝑊 − 𝑀 𝑀𝑇 𝑀 −1 𝑀𝑇 𝑊 2
𝑀
The explicit expression will be very complicated!
We employ Black-box optimization of the function,𝑓 𝑀 = 𝑊 − 𝑀 𝑀𝑇 𝑀 −1 𝑀𝑇 𝑊 2(Numerical)
Experiment

data acquisition

data modeling𝑧𝑖 = 𝑓 𝑀𝑖

𝑧=𝑔 𝑀

candidate generation

MachineLearning

QUBO form

𝑀𝑖+1 = argmin𝑀 𝑔(𝑀)Optimization
量子アニーリングと機械学習の交差点 / 2021/11/12 @ IBIS2021© DENSO CORPORATION All Rights Reserved.

17 / 39min 𝑓 𝑀𝑖Results (8 x 100 => 8 x 3 & 3 x 100 [N=24])𝑖≤𝑛gBOCS: tuned BOCS-SAvBOCS: vanilla BOCS-SA
FMQA08: FMQA (k=8)
FMQA12: FMQA (k=12)

Original repeated
1-rank approximation
Ground state(Exact solution)𝑛 (# of acquired data)

量子アニーリングと機械学習の交差点 / 2021/11/12 @ IBIS2021© DENSO CORPORATION All Rights Reserved.

18 / 39

Quality of solutions and execution time
10 randomly generated instances
25 independent runs

random gBOCS vBOCS

FMQA08 FMQA12GS count36265853

GS win count05042

Energy win count070210.355124394

3719

3714

Time

1/480
1/70

Comparable/better solutions with faster execution time
量子アニーリングと機械学習の交差点 / 2021/11/12 @ IBIS2021© DENSO CORPORATION All Rights Reserved.

19 / 39

BOCS-SA and BOCS-QA

N=64𝑖≤𝑛min 𝑓 𝑀𝑖𝑖≤𝑛min 𝑓 𝑀𝑖

N=24𝑛 (# of acquired data)𝑛 (# of acquired data)

We see some differences between SA and QA, but not significant
量子アニーリングと機械学習の交差点 / 2021/11/12 @ IBIS2021© DENSO CORPORATION All Rights Reserved.

20 / 39

Vibration-resistant design (VRD)
Mounting holes

Power control unit
https://www.denso.com/jp/ja/news/newsroom/2019/20190522-01//media/8c44d0012d0148f893903eb515d52b67.ashx

Printed circuit board (TOYOTA Yaris)

Vibrations in the power control unit will damage the printed circuit board.
To reduce the damage, we design a PCB with a high resonance frequency.
That generally causes many mounting holes. The goal is to design a PCB
with a higher resonance frequency and fewer mounting holes.

量子アニーリングと機械学習の交差点 / 2021/11/12 @ IBIS2021© DENSO CORPORATION All Rights Reserved.

21 / 39

Formulation of the problem

Collaboration with T. Matsumori and M. Taki

Candidate points (𝑥𝑖 = 0)
Selected points (𝑥𝑖 = 1)
Keep out area
Minimize # of mounting holes
Resonance frequency 𝒢 𝑥 constraint
by finite element method (FEM)
Modified cost function:
Higher resonance frequency
Target # of mounting holes

Non-linear & non-QUBO combinatorial optimization problem
量子アニーリングと機械学習の交差点 / 2021/11/12 @ IBIS2021© DENSO CORPORATION All Rights Reserved.

22 / 39

How to solve it?

The problem to be solved:
argmin 𝓌𝐹 𝐹 𝒙 − 𝑃ത 2 − 𝓌𝐺 𝒢 𝒙𝒙∈ 0,1 𝑛𝑓 𝒙 = 𝓌𝐹 𝐹 𝒙 − 𝑃ത 2 − 𝓌𝐺 𝒢 𝒙

𝑓 𝒙𝑖(Numerical)
Experiment

data acquisition

data modeling𝑧𝑖 = 𝑓 𝒙𝑖

𝑧=𝑔 𝒙

candidate generation

𝒙𝑖+1

QUBO form

𝒙𝑖+1 = argmin𝒙 𝑔(𝒙)Optimization

量子アニーリングと機械学習の交差点 / 2021/11/12 @ IBIS2021© DENSO CORPORATION All Rights Reserved.

MachineLearning

23 / 39

Results (N=12) [𝟐𝟏𝟐 = 𝟒𝟎𝟗𝟔]

SA 3p: optimal + 2 neighbors
SA 1p: optimal

Solution quality
Model accuracy

SA and QA outperform random search (RS).

量子アニーリングと機械学習の交差点 / 2021/11/12 @ IBIS2021© DENSO CORPORATION All Rights Reserved.

24 / 39

1. ブラックボックス最適化••

非可逆行列圧縮
基板設計最適化

2. 量子アニーリング•

次世代量子アニーラ向けアルゴリズム

量子アニーリングと機械学習の交差点 / 2021/11/12 @ IBIS2021© DENSO CORPORATION All Rights Reserved.

25 / 39

Solving combinatorial optimization problems
by qubits• Quantum Annealing (QA)
• Adiabatic Quantum Computation (AQC)• Quantum Approximate Optimization Algorithm (QAOA)• Variational Quantum Eigensolver (VQE)
Diabatic vs. Adiabatic, Global vs. Local search,
Unitary transformations (types & control parameters), etc…
Shorter processing times are desired in NISQ devices. We focus on dynamicalaspects of QA with non-stoquastic Hamiltonian. A counter-diabatic driving andvariational optimization are studied. We propose a new algorithm forcombinatorial optimization. Comparison with conventional QA and simulated
annealing (SA) will be presented, as well as a proposed oracle is discussed.

量子アニーリングと機械学習の交差点 / 2021/11/12 @ IBIS2021© DENSO CORPORATION All Rights Reserved.

26 / 39

Model Hamiltonian
with counter-diabatic drivingℋ = 𝐴 𝑡 ℋ𝑧 + 𝐵 𝑡 ℋ𝑥 + ℋ𝑦

断熱

非断熱

1σ𝑖<𝑗 𝜎𝑖𝑧 𝜎𝑗𝑧𝑁−1ℋ𝑥 = − σ𝑖 𝜎𝑖𝑥: uniform transverse field𝑦σℋ𝑦 = − 𝑖 𝐶𝑖 𝑡 𝜎𝑖: non-uniform y-fieldℋ𝑧 = −

where the coefficients are

CD制御: cost function(non-stoquastic, counter-diabatic)

𝐴 𝑡 = 𝑎𝑡/𝜏𝑎, 𝑏, 𝑐𝑖 : parameters

𝐵 𝑡 = 𝑏 1 − 𝑡/𝜏𝜏 : annealing time

𝐶𝑖 𝑡 = 𝑐𝑖 sin2 𝜋𝑡/𝜏(𝜏 = 𝑎 = 1 for simplicity)
Sels 2017, Takahashi 2017, Prielinger 2021

量子アニーリングと機械学習の交差点 / 2021/11/12 @ IBIS2021© DENSO CORPORATION All Rights Reserved.

27 / 39

Time evolution of the system𝑑ȁ𝜓ۧ = −𝑖ℋ ȁ𝜓ۧ = −𝑖 𝐴 𝑡 ℋ𝑧 + 𝐵 𝑡 ℋ𝑥 + ℋ𝑦 ȁ𝜓ۧ𝑑𝑡
Initial : transverse field (makes superposition of all possible states)
Intermediate :(1) transverse field is decreasing,(2) effective local field* in z is growing,(3) non-uniform y-field is applied (we assume it’s dominant.)
Final : effective local field
Initialȁ0ۧIntermediate

Finalห𝜓Final ൿห𝜓Initial ൿyx
ȁ1ۧ* This caused by spin-spin interactions in the cost function
量子アニーリングと機械学習の交差点 / 2021/11/12 @ IBIS2021© DENSO CORPORATION All Rights Reserved.

28 / 39

Mean-field model & variational optimization
Mean field approximation:ℋ = −𝑡 𝜎 𝑧 𝜎 𝑧 − 𝑏 1 − 𝑡 𝜎 𝑥 − 𝑐 sin2 𝜋𝑡 𝜎 𝑦𝜎𝑧 = 𝜓 𝜎𝑧 𝜓optimize magnetization 𝜎 𝑧ℋ𝑥𝑏𝑐𝑏𝑐𝑏opt = 0.539
𝑐opt = 1.565(𝜋Τ2 = 1.571)

𝑐𝑏t

量子アニーリングと機械学習の交差点 / 2021/11/12 @ IBIS2021© DENSO CORPORATION All Rights Reserved.

29 / 39

Measures of optimizationin many-body systems
We use two measures:• Fidelity(𝑃GS ) : probability to find the ground state from the final state2𝑃GS = 𝜓Final 𝜓GS , which can’t be observed by experiment.
A oracle to measure this probability is required.
(This requirement will be relaxed in later analysis.)
Error = 1 − 𝑃GS• Energy : energy of the final state 𝜓Final ℋ𝑧 𝜓Final , which canbe observed.
Ground states and excited states are not distinguishable.

量子アニーリングと機械学習の交差点 / 2021/11/12 @ IBIS2021© DENSO CORPORATION All Rights Reserved.

30 / 39

Variational optimization(Ferromagnetic, fully connected, 𝑵 = 𝟐𝟎)
fidelityand energy.

Robust:• transverse fieldSensitive: (reciprocal)• non-uniform field• annealing time★
(.54, 1.56)

annealing timeoptimized by fidelity

transverse field

Parameters can beEnergy

non-uniform field
量子アニーリングと機械学習の交差点 / 2021/11/12 @ IBIS2021© DENSO CORPORATION All Rights Reserved.

31 / 39

Variational optimization (size dependence)
Ferromagnetic, 𝑵 = 𝟐~𝟐𝟎
𝑁𝑏opt

𝑁𝑐opt
𝑁𝑁
use 𝑏opt& 𝑐optforfurther analyses

N

N
Error (= 1 − 𝑃GS)

size

Scaling relation
between size andoptimization error
is observed.

t
量子アニーリングと機械学習の交差点 / 2021/11/12 @ IBIS2021© DENSO CORPORATION All Rights Reserved.

N32 / 39Sequential optimization (Ferromagnetic)Now, we need a little relaxed
oracle to answer gradients(+/-) at𝑐𝑖 = 0.

量子アニーリングと機械学習の交差点 / 2021/11/12 @ IBIS2021© DENSO CORPORATION All Rights Reserved.

Sequential estimation for N=8# of fixed parameters

In general systems (e.g., randomsystems), we need to identify eachparameter’s sign (+ or −). For thispurpose, one point estimation of
the gradient at 𝑐𝑖 = 0 is enough toidentify the sign.
signs of gradients= sign of parameters= spin orientations in GS

33 / 39Sequential QGOparameters•transverse field•non-uniform field in Y

sensitivity analysis at 𝑐𝑖 = 0

量子アニーリングと機械学習の交差点 / 2021/11/12 @ IBIS2021© DENSO CORPORATION All Rights Reserved.

34 / 39Sequential QGO of random systemsℋ = 𝐴 𝑡 ℋ𝑧 + 𝐵 𝑡 ℋ𝑥 + ℋ𝑦
ℋ𝑧 = − σ𝑖<𝑗 𝐽𝑖𝑗 𝜎𝑖𝑧 𝜎𝑗𝑧: cost functionℋ𝑥 = − σ𝑖 𝜎𝑖𝑥: uniform transverse field𝑦: non-uniform y-fieldℋ𝑦 = − σ𝑖 𝐶𝑖 𝑡 𝜎𝑖This model is equivalent to Max-cut problem.
(𝐽𝑖𝑗 = −𝑤𝑖𝑗 ; 𝑤𝑖𝑗 is the weight matrix of the problem)

We employ SK model interaction:
𝑃 𝐽𝑖𝑗 ∼12𝜋𝜎exp −22𝐽𝑖𝑗

2𝜎

2 =,where𝜎2

量子アニーリングと機械学習の交差点 / 2021/11/12 @ IBIS2021© DENSO CORPORATION All Rights Reserved.

1𝑁−135 / 39

Sequential QGO of random systems
Spin index #-fidelity
optimization

Success ratio
100%# of fixed parameters++testingfixingfixed+-+---+-+-

non-uniform y-field
Fidelity optimization always finds the ground state.
量子アニーリングと機械学習の交差点 / 2021/11/12 @ IBIS2021© DENSO CORPORATION All Rights Reserved.

36 / 39Sequential QGO of random systems
Spin index #-Energyoptimization

Success ratio??%# of fixed parameters-testing
fixingfixed+-++---+-+-

non-uniform y-field
Energy optimization also finds the ground state.
量子アニーリングと機械学習の交差点 / 2021/11/12 @ IBIS2021© DENSO CORPORATION All Rights Reserved.

37 / 39Comparison with SA and QA
Short time (anneal time 𝜏 = 1, 5) QA does not find the exact solution, while SA andQGO have a certain volume of basin.

N
100 random instances are analyzed and their performances is averaged.
While sequential QGO runs N(N+3)/2 QA tasks, no adjustment is applied.

QGO employs non-uniform y-field to probe the ground state or low energystates enriched during annealing process.

量子アニーリングと機械学習の交差点 / 2021/11/12 @ IBIS2021© DENSO CORPORATION All Rights Reserved.

38 / 39

0. 量子アニーリング？
1. ブラックボックス最適化••

非可逆行列圧縮
基板設計最適化

2. 量子アニーリング•

機械学習 > 量子 な話

量子 > 機械学習 な話

次世代量子アニーラ向けアルゴリズム

皆さん一緒に、離散変数のブラックボックス最適化やりませんか？

量子アニーリングと機械学習の交差点 / 2021/11/12 @ IBIS2021© DENSO CORPORATION All Rights Reserved.

39 / 39