ラック・セキュリティごった煮ブログ

セキュリティエンジニアがエンジニアの方に向けて、 セキュリティやIT技術に関する情報を発信していくアカウントです。

【お知らせ】2021年5月10日~リニューアルオープン!今後はこちらで新しい記事を公開します。

株式会社ラックのセキュリティエンジニアが、 エンジニアの方向けにセキュリティやIT技術に関する情報を発信するブログです。(編集:株式会社ラック・デジタルペンテスト部)
当ウェブサイトをご利用の際には、こちらの「サイトのご利用条件」をご確認ください。

デジタルペンテスト部提供サービス:ペネトレーションテスト

SECCON CTF 13国内決勝 Writeup

デジタルペンテスト部の中嶋です。 2025年の3月1日(土)から3月2日(日)にかけて開催されたSECCON CTF 13の国内決勝にチームBunkyoWesterns、プレイヤー名parabola0149として参加しました。 実はSECCON CTF 13の予選では別のチームとして参加して予選を通過することができなかったのですが、BunkyoWesternsのメンバーから声がかかり国内決勝に出場できることになりました。 私が問題に取り組んでいる間にチームのメンバーがモリモリ得点を稼いでいて、国内決勝で優勝することができました。

SECCONの国内決勝に参加する機会を与えてくれて、共に国内決勝に参加したBunkyoWesternsのみなさん、ありがとうございました。

SECCON CTF 13国内決勝の概要

SECCON CTF 13の国内決勝は2日間開催されていて、日中の競技時間にCTFのサーバーにアクセスすることができます。 出題された問題の形式としてはJeopardy形式とKing of the Hill形式の2種類がありました。

Jeopardyは通常のCTFの問題形式で、問題を解いてフラグを提出すると得点が入る形式です。 Jeopardy形式の問題は配布ファイルをダウンロードすれば夜間にも問題に取り組めるようになっていました。 King of the Hillの詳細はHECCONという問題の解説で説明しますが、問題サーバーにアクセスする必要があるので日中の競技時間にしかできず、さらに競技時間中は付きっきりで取り組む必要のある形式でした。 そこで、CTFの競技時間はKing of the Hill形式の問題に取り組んで夜間にJeopardy形式の問題に取り組む方針としました。

私はCryptoが得意ジャンルであり、King of the Hill形式の問題の中でHECCONという問題がCryptoっぽい問題だったので日中はその問題に取り組んでいました。 HECCONのスコアの合計や順位はスコアサーバ等には載っていなかったと思うのですが、手元で計算して確認したところスコアの合計が2122 ptで、9チーム中2位になることができました。

夜間はJeopardy形式のCryptoの問題であるx_xを解きました。 国内決勝に参加したチームの中ではx_xを解けたのは私たちのチームだけでした。 Jeopardy形式のCryptoの問題は他にもあり特にhell_summonという問題は私も取り組んでいたのですが、なかなか解けないでいたところチームのメンバーが解いてくれたので助かりました。

私が解いた問題はx_xHECCONの2問なので、この2問についてWriteupで解説します。

SECCON CTF 13決勝開催中に何枚か写真を取っていたので紹介します。

今回のCTFではチームBunkyoWesternsとして参加することになりました

SECCON CTF 13の決勝ではハードウェア問も出題されました

SECCON CTF 13決勝のチーム一覧です。国内決勝と国際決勝が同じ会場で行われていました

kurenaifさんとアルパカのパネルがありました。同じ会場でAlpacaHack開催の決勝観戦CTFがあったみたいです

フラグを提出するとスクリーンにkurenaifさんが現れます

[crypto] x_x (1 Solved, 500 pt)

x_x < hello

nc x-x.dom.seccon.games 31180

x_x.tar.gz 274320b9a95ee5eaa878330d04bbf4eaeeb9da76

問題の概要

この問題では以下のファイルが配布されています。

  • 問題サーバーで動かしているSageMathプログラム (server.sage)
  • 問題サーバーを動かすときに使われているDockerfile (Dockerfile)

この問題では q = 31, n = 180, v = 165, o = n - v がパラメータとして使用されています。また、特に断りがなければこの問題で使用する数は位数が q = 31 の有限体 \mathbb{F} _ {q} の要素とします。

問題サーバーに接続すると最初にランダムな n \times n 行列 S  0 \le i \lt o としてランダムな v \times o 行列 A _ i 、ランダムな v \times v 行列 B _ i が生成されます。  O を零行列としたときに、以下の n \times n 行列を F _ i としています。

 \begin{align}
F _ i &= \begin{bmatrix}
O & A _ i^ \mathsf{T} \\
A _ i & B _ i
\end{bmatrix}
\end{align}

 n \times n 行列 P _ i  P _ i = S^ \mathsf{T} F _ i S として計算されて、公開鍵として表示されます。

b"give me flag"というデータをhash関数に入力して h _ i (問題のプログラムではdata)が計算されています。  o 次元のベクトル s (問題のプログラムではsig)を入力して、 s P _ i s^ \mathsf{T} = h _ i が成り立つようにするとフラグが表示されます。

  • server.sage
from Crypto.Util.number import bytes_to_long
from Crypto.Cipher import AES
import secrets
import os
import signal
import hashlib

signal.alarm(300)

flag = os.getenv('flag', "SECCON{sample}")

q=31
n=180
v=165
o=n-v
K = GF(q)

def gen_key():
    S = random_matrix(K, n, n)
    # mats = matrix.identity(K, n)

    Bs = [random_matrix(K, v, v)]
    As = [random_matrix(K, v, o)]

    P = [[0 for _ in range(o)] for _ in range(o)]
    for i in range(o):
        P[i][i] = 1
    P = P[1:] + P[:1]
    P = matrix(K, P)

    for i in range(o-1): 
        Bs.append(random_matrix(K, v, v))
        As.append(As[-1] * P)

    Fs = []
    for i in range(o):
        F = block_matrix(K, [
            [0, (As[i]).transpose()],
            [As[i], Bs[i]]
        ])
        Fs.append(F)

    Ps = []
    for F in Fs:
        Ps.append(S.transpose() * F * S)
    return (Ps, (S, Fs))

def verify(pubkey, sig, data):
    return [sig * P * sig for P in pubkey] == data

def hash(data):
    res = []
    while len(res) < o:
        data = hashlib.sha512(data).digest()
        val = bytes_to_long(data)
        while val > 0:
            res.append(val % q)
            val = int(val / q)
            if len(res) == o:
                break
    assert len(res) == o
    return res

pubkey, _privkey = gen_key()
print(f"pubkey={[list(mat) for mat in pubkey]}")
sig = input("sig=")[1:-1] #format: "(x0,x1,x2,x3, ... ,xo)"
sig = vector(map(K, sig.split(",")))
data = hash(b"give me flag")
if verify(pubkey, sig, data):
    print(flag)
  • Dockerfile
FROM sagemath/sagemath:10.5

USER root
COPY server.sage .
ENV flag="SECCON{sample}"
RUN sage --pip install pycryptodome
RUN apt-get update && apt-get install socat -y
CMD ["sage", "solver.sage"]
CMD ["socat", "TCP-L:31180,fork,reuseaddr", "EXEC:'sage server.sage'"]

解法

この問題では s P _ i s^ \mathsf{T} = h _ i という形式の方程式が o 個あるので、 o 個の変数 x _ i を用意すれば全ての方程式を満たすようなベクトル s が求まりそうです。 最初に試したのは、 s の成分の内の o 個をそれぞれ変数 x _ i として残りを 0 としたときの式を計算して、グレブナー基底を使用してそれらの式を簡約化して解くことができないかと考えました。

 \begin{align}
s &= \left[ \begin{array}{cccc|cccc}
x _ 0 & x _ 1 & \cdots & x _ {o - 1} & 0 & 0 & \cdots & 0
\end{array} \right]
\end{align}

この場合、 s P _ i s^ \mathsf{T} = h _ i は2次の方程式になります。 実際にこの方法を試したところグレブナー基底の計算に時間がかかって完了せず、さすがにこの方法では解くことができなそうでした。

  • groebner_basis.sage
#!/usr/bin/env sage

from pwn import remote
from Crypto.Util.number import bytes_to_long
import hashlib
from ast import literal_eval


host = 'x-x.dom.seccon.games'
# host = 'localhost'
port = 31180


q=31
n=180
v=165
o=n-v
K = GF(q)

def hash(data):
    res = []
    while len(res) < o:
        data = hashlib.sha512(data).digest()
        val = bytes_to_long(data)
        while val > 0:
            res.append(val % q)
            val = int(val / q)
            if len(res) == o:
                break
    assert len(res) == o
    return res


with remote(host, port) as conn:
    data = conn.recvline()
    pubkey = [matrix(K, i) for i in literal_eval(data.splitlines()[0].decode().strip('pubkey='))]
    conn.recvuntil(b'sig=')
    
    h = hash(b"give me flag")
    P = PolynomialRing(K, ','.join(f'x{i}' for i in range(o)))
    sig_ = vector(list(P.gens()) + [0] * v)
    f_list = [sig_ * P * sig_ - y for P, y in zip(pubkey, h)]
    I = Ideal(f_list)
    B = I.groebner_basis()
    print([i.monomials() for i in B])

問題を確認すると行列 F _ i には o \times o の零行列が含まれていて、これを上手く使えば s P _ i s^ \mathsf{T} = h _ i を1次の方程式にして簡単に解くことができるようになると考えました。 ただ、行列 P _ i  P _ i = S^ \mathsf{T} F _ i S というように行列 S が掛けられていて零行列の部分が無くなっているので、一旦 S が無い場合の問題サーバーを再現してこの方法を試すことにしました。

先ほどと同様に s の成分の内の o 個をそれぞれ変数 x _ i としました。 ただし、残りを 0 とすると s P s^ \mathsf{T} = 0 となってしまうので、代わりに残りの成分を 1 としました。

 \begin{align}
s &= \left[ \begin{array}{cccc|cccc}
x _ 0 & x _ 1 & \cdots & x _ {o - 1} & 1 & 1 & \cdots & 1
\end{array} \right]
\end{align}

 S が無い場合はこの方法で条件を満たす s を求めることができました。 以下のプログラムではグレブナー基底を使用したときの名残でI = Ideal(f_list); V = I.variety()として連立方程式の解を求めていますが、1次の連立方程式なので普通の行列演算で解を求めることができるはずです。

  • server.sage の変更点
$ diff server.sage server_2.sage
45c45
<         Ps.append(S.transpose() * F * S)
---
>         Ps.append(F)
  • without_S_matrix.sage
#!/usr/bin/env sage

from pwn import remote
from Crypto.Util.number import bytes_to_long
import hashlib
from ast import literal_eval


host = 'x-x.dom.seccon.games'
# host = 'localhost'
port = 31180


q=31
n=180
v=165
o=n-v
K = GF(q)

def hash(data):
    res = []
    while len(res) < o:
        data = hashlib.sha512(data).digest()
        val = bytes_to_long(data)
        while val > 0:
            res.append(val % q)
            val = int(val / q)
            if len(res) == o:
                break
    assert len(res) == o
    return res


with remote(host, port) as conn:
    data = conn.recvline()
    pubkey = [matrix(K, i) for i in literal_eval(data.splitlines()[0].decode().strip('pubkey='))]
    conn.recvuntil(b'sig=')
    
    h = hash(b"give me flag")
    P = PolynomialRing(K, ','.join(f'x{i}' for i in range(o)))
    sig_ = vector(list(P.gens()) + [1] * v)
    f_list = [sig_ * P * sig_ - y for P, y in zip(pubkey, h)]
    I = Ideal(f_list)
    V = I.variety()
    sig = vector([V[0][var] for var in P.gens()] + [1] * v)
    
    conn.sendline(str(sig).encode())
    data = conn.recvall()
    print(data)

この方針で問題が解けそうなので、次は S を打ち消す方法や S が無い場合のように o \times o の零行列の箇所を作る方法を考えました。 もう少し言うと、 C P _ i C^ \mathsf{T} = O が成り立つような o \times n 行列 C を求めたいです。

 P _ i の代わりに F _ i を使って考えると、任意の o \times o 行列 X と零行列 O を連結した行列を C = \left[ \begin{array}{c}
X
\end{array} \:\:\: \begin{array}{c}
O
\end{array} \right] とすると、

 \begin{align}
C F _ i C^ \mathsf{T} &= \begin{bmatrix}
X & O
\end{bmatrix} \begin{bmatrix}
O & A _ i^ \mathsf{T} \\
A _ i & B _ i
\end{bmatrix} \begin{bmatrix}
X & O
\end{bmatrix}^ \mathsf{T} \\
&= \begin{bmatrix}
X & O
\end{bmatrix} \begin{bmatrix}
O & A _ i^ \mathsf{T} \\
A _ i & B _ i
\end{bmatrix} \begin{bmatrix}
X \\
O
\end{bmatrix} \\
&= \begin{bmatrix}
O & X A _ i^ \mathsf{T}
\end{bmatrix} \begin{bmatrix}
X \\
O
\end{bmatrix} \\
&= O
\end{align}

となるので C F _ i C^ \mathsf{T} = O が成り立ちます。 このような C  F _ i の具体的な成分の情報から推測せずに計算で求める方法を考えます。

行列 C  o 列目以降の成分が全て 0 なので、 F _ i と掛けると B _ i の情報が無くなります。 そこで、 F _ i から B _ i の情報のみを取り出して、それに対して掛けると 0 になるベクトルの集合(カーネル)を求めると C が得られるのではないかと考えました。  B _ i の情報を取り出すために、 F _ i - F _ i^ \mathsf{T} という式を考えると以下のようになります。

 \begin{align}
F _ i - F _ i^ \mathsf{T} &= \begin{bmatrix}
O & A _ i^ \mathsf{T} \\
A _ i & B _ i
\end{bmatrix} - \begin{bmatrix}
O & A _ i^ \mathsf{T} \\
A _ i & B _ i
\end{bmatrix}^ \mathsf{T} \\
&= \begin{bmatrix}
O & A _ i^ \mathsf{T} \\
A _ i & B _ i
\end{bmatrix} - \begin{bmatrix}
O & A _ i^ \mathsf{T} \\
A _ i & B _ i^ \mathsf{T}
\end{bmatrix} \\
&= \begin{bmatrix}
O & O \\
O & B _ i - B _ i^ \mathsf{T}
\end{bmatrix}
\end{align}

 C  F _ i - F _ i^ \mathsf{T} に掛けると、

 \begin{align}
C \left( F _ i - F _ i^ \mathsf{T} \right) &= \begin{bmatrix}
X & O
\end{bmatrix} \begin{bmatrix}
O & O \\
O & B _ i - B _ i^ \mathsf{T}
\end{bmatrix} \\
&= \begin{bmatrix}
O & O
\end{bmatrix}
\end{align}

というように零行列になるので、 C から生成されるベクトルは全て F _ i - F _ i^ \mathsf{T} カーネルに含まれることが分かります。  F _ i - F _ i^ \mathsf{T} カーネルを確認すると C に関係ないものも含まれていますが、 0 \le i \lt o に対する F _ i - F _ i^ \mathsf{T} カーネルを計算して全てに共通する部分空間を求めると C を求めることができました。

ここまで F _ i を使って考えてきましたが、 S 逆行列が存在する場合は P _ i に対しても同じ計算で C P _ i C^ \mathsf{T} = O を満たすような o \times n 行列 C を求めることができます。

 C を求めた後は、ランダムな v \times n 行列 R を生成して C, R を連結した行列 D = \begin{bmatrix}
C \\
R
\end{bmatrix} を用意します。  D P _ i D^ \mathsf{T} を計算すると、

 \begin{align}
D P _ i D^ \mathsf{T} &= \begin{bmatrix}
C \\
R
\end{bmatrix} P _ i \begin{bmatrix}
C \\
R
\end{bmatrix}^ \mathsf{T} \\
&= \begin{bmatrix}
C \\
R
\end{bmatrix} P _ i \begin{bmatrix}
C^ \mathsf{T} & R^ \mathsf{T}
\end{bmatrix} \\
&= \begin{bmatrix}
C P _ i \\
R P _ i
\end{bmatrix} \begin{bmatrix}
C^ \mathsf{T} & R^ \mathsf{T}
\end{bmatrix} \\
&= \begin{bmatrix}
C P _ i C^ \mathsf{T} & C P _ i R^ \mathsf{T} \\
R P _ i C^ \mathsf{T} & R P _ i R^ \mathsf{T}
\end{bmatrix} \\
&= \begin{bmatrix}
O & C P _ i R^ \mathsf{T} \\
R P _ i C^ \mathsf{T} & R P _ i R^ \mathsf{T}
\end{bmatrix}
\end{align}

となるので、 F _ i のように o \times o の零行列が含まれています。 よって、 S が無い場合と同じ方法で s \left( D P _ i D^ \mathsf{T} \right) s^ \mathsf{T} = h _ i を満たす s を求めることができます。  s D を問題サーバーに入力すると、

 \begin{align}
\left( s D \right) P _ i \left( s D \right)^ \mathsf{T} &= s D P _ i D^ \mathsf{T} s^ \mathsf{T} \\
&= s \left( D P _ i D^ \mathsf{T} \right) s^ \mathsf{T} \\
&= h _ i
\end{align}

となるので、問題の条件を満たしてフラグを取得することができます。

  • solve.sage
#!/usr/bin/env sage

from pwn import remote
from Crypto.Util.number import bytes_to_long
import hashlib
from ast import literal_eval
import functools


host = 'x-x.dom.seccon.games'
# host = 'localhost'
port = 31180


q=31
n=180
v=165
o=n-v
K = GF(q)

def hash(data):
    res = []
    while len(res) < o:
        data = hashlib.sha512(data).digest()
        val = bytes_to_long(data)
        while val > 0:
            res.append(val % q)
            val = int(val / q)
            if len(res) == o:
                break
    assert len(res) == o
    return res


with remote(host, port) as conn:
    data = conn.recvline()
    pubkey = [matrix(K, i) for i in literal_eval(data.splitlines()[0].decode().strip('pubkey='))]
    conn.recvuntil(b'sig=')
    
    h = hash(b"give me flag")
    C = functools.reduce(lambda x, y: x & y, [(i - i.T).kernel() for i in pubkey]).matrix()
    assert all((C * i * C.T).is_zero() for i in pubkey)
    P = PolynomialRing(K, ','.join(f'x{i}' for i in range(o)))
    while True:
        D = block_matrix(K, [[C], [random_matrix(K, v, n)]])
        sig_ = vector(list(P.gens()) + [1] * v) * D
        f_list = [sig_ * P * sig_ - y for P, y in zip(pubkey, h)]
        I = Ideal(f_list)
        V = I.variety()
        if V:
            break
    sig = vector([V[0][var] for var in P.gens()] + [1] * v) * D
    
    conn.sendline(str(sig).encode())
    data = conn.recvall()
    print(data)
SECCON{x_x=x.x=x~x=xxx=x_x}

[King of the Hill] HECCON (2122 pt, Rank: 2)

Homomorphic Encryption Coding CONtest

First download file (including pubkey and some useful information etc):

https://static.final.seccon.jp/heccon.tar.gz

then go to the following website:

https://heccon.dom.seccon.games

問題の概要

この問題はKing of the Hill形式の問題で、通常のJeopardy形式の問題とはタイプが異なります。 問題文に書かれたURLにアクセスすると、以下のページが表示されました。

HECCONの問題の概要ですが、この問題は暗号化したまま加算や乗算ができる準同型暗号についての問題です。 HECCONで使用されている暗号の特徴は以下の通りです。

  • 暗号化したまま加算と乗算が両方とも可能
  • 乗算の回数には制限がある
  • 除算や平方根などの他の計算は(簡単には)できない

この問題では最初にランダムな値 x が生成されます。  x に対してある計算がなされて、その出力を y = f \! \left( x \right) としています。 入力 x が暗号化されて配布されるので、暗号化された状態で可能な操作のみを使用して y' = f' \! \left( x \right) を計算します。  y' を問題サーバに送信すると、 y' が復号されて y  y' が比較されます。  y  y' の平均絶対誤差(Mean Absolute Error; MAE)が小さいほど順位が高くなります。 つまり問題サーバで何らかの計算が行われるので、暗号化したままなるべく似たような計算をしたいというのがこの問題でやりたいことです。

HECCONはKing of the Hill形式の問題であり、そのルールをざっくりと説明すると、1ラウンド5分で各ラウンドごとに順位に応じたスコアが加算されます。 つまりこのルールでは高い順位を取り続けることが重要になるので、King of the Hill形式の問題は競技時間中ずっと取り組み続ける必要がありました。 HECCONでは平均絶対誤差が小さいほど順位が高くなり多くのスコアを得ることができます。

また、HECCONには42ラウンドのステージが4つあります。 各ステージごとに出力 y = f \! \left( x \right) の計算式が異なり、ステージごとにどのように計算するかを考えなければならないので、実質4問分の問題を解いているような感覚でした。

この問題では以下のファイルが配布されています。

  • 問題と配布ファイルの説明 (README.md)
  • 公開鍵と秘密鍵を生成するPythonファイル (generate_key.py)
  • 公開鍵 (pubkey)
  • 問題を解くときに使えるDockerfile (Dockerfile)
  • 問題を解くときに使えるコードが書かれたファイル*1 (snippet.md)

さらに各ステージごとに以下のファイルが配布されます。

  • 入力 x の生成や出力 y = f \! \left( x \right) の計算、 y  y' の平均絶対誤差の計算を行うプログラム (challenge.py)
  • 暗号化された入力 x (enc)

ステージ共通の配布ファイルの内容はそれぞれ以下の通りです。

  • README.md
# HECCON dist files

This directory contains the following files:

- Dockerfile: The docker environtment you can use for this challenge
- pubkey: The public key for homomorphic encryption
- generate_key.py: The script to generate the public and private keys (only the public key is distributed as pubkey)
- snippet.md: Some useful snippets to for this challenge

Note that `helayers` doesn't work with recent versions of Python. I confirmed that it works with Python 3.11.

Don't forget to download challenge files from the challenge website.
  • generate_key.py
import pyhelayers

pubkey_path = "/path/to/pubkey"
privkey_path = "/path/to/privkey"

req = pyhelayers.HeConfigRequirement(
    num_slots=2**13,
    integer_part_precision=10,
    fractional_part_precision=30,
    multiplication_depth=10,
    security_level=128,
)
he_context = pyhelayers.SealCkksContext()
he_context.init(req)
he_context.save_to_file(pubkey_path)
he_context.save_secret_key_to_file(privkey_path)
  • Dockerfile
FROM python:3.11.11-slim

RUN pip3 install --no-cache-dir pyhelayers==1.5.4.0 numpy==2.2.2
  • snippet.md
# Snippet

You can use the following codes as a snippet.

```python
from pathlib import Path
import numpy as np
import pyhelayers

pubkey_path = "/path/to/pubkey"
privkey_path = "/path/to/privkey"
```

Generation of your own pubkey-privkey pair for test.
I recommend to generate it because you can debug your code easily.

```python
req = pyhelayers.HeConfigRequirement(
    num_slots=N,
    integer_part_precision=10,
    fractional_part_precision=30,
    multiplication_depth=10,
    security_level=128,
)
he_context = pyhelayers.SealCkksContext()
he_context.init(req)
he_context.save_to_file(pubkey_path)
he_context.save_secret_key_to_file(privkey_path)
```

Loading of key

```python
he_context = pyhelayers.SealCkksContext()
he_context.load_from_file(pubkey_path)
# If you have privkey
# he_context.load_secret_key_from_file(privkey)
```

Creation of encoder

```python
encoder = pyhelayers.Encoder(he_context)
```

Encryption

```python
x = np.linspace(-1, 1, 2**13)
cx = encoder.encode_encrypt(x)
# You can also encrypt empty array
cy = encoder.encode_encrypt([])
```

Decryption

```python
x = encoder.decrypt_decode_double(cx)
```

Loading of encryption data from file

```python
cx = encoder.encode_encrypt([])
cx.load_from_file("/path/to/enc")
```

Save of encryption data to file
This type of files can be submitted and the server evaluate the score.

```python
cx.save_to_file("/path/to/submission")
```

Addition
Note that the operation is in-place.

```python
# vector
cx.add(cy)
# scalar
cx.add_scalar(1.0)
```

Subtraction

```python
# vector
cx.sub(cy)
```

Multiplication
Note that you can't multiply vectors so much. See `multiplication_depth` in HeConfigRequirement.

```python
# vector
cx.multiply(cy)
# scalar
cx.multiply_scalar(2.0)
```

Check current chain index

```python
cx.get_chain_index()
```

Some useful operations
Ref: https://ibm.github.io/helayers/reference/gen_pages/FunctionEvaluator.html

```python
fe = pyhelayers.FunctionEvaluator(he_context)
```

Stage 1 (Score: 425, Rank: 5)

ここから各ステージの説明ですが、最初に順位表があった方が各ステージの展開が分かりやすい気がしたので先に載せておきます。 順位表の画像は運営の方がXに載せていたものです(https://x.com/y011d4/status/1896198494760939795)。 グラフの縦軸が平均絶対誤差なので、下に行くほど順位が高くなります。

ステージ1のchallenge.pyは以下の通りです。

  • challenge.py
from pathlib import Path

import numpy as np
import pyhelayers


N = 2**13


class Challenge:
    @classmethod
    def generate(cls, pubkey: Path, dist_dir: Path):
        """Generate the challenge data.

        dist_dir is a distributed directory to contestants.

        Args:
            pubkey (Path): Path to the public key file.
            dist_dir (Path): Path to the directory where the challenge data is saved.
        """
        he_context = pyhelayers.SealCkksContext()
        he_context.load_from_file(str(pubkey))
        encoder = pyhelayers.Encoder(he_context)
        x = np.random.uniform(-3, 3, N)
        cx = encoder.encode_encrypt(x)
        cx.save_to_file(str(dist_dir / "enc"))
        (dist_dir / "challenge.py").write_text(Path(__file__).read_text())

    @classmethod
    def evaluate(
        cls, pubkey: Path, privkey: Path, enc: Path, submission: Path
    ) -> float:
        """Evaluate the submission.

        submission is a path to a submitted file and the server evaluate it.
        Contestants are supposed to get better MAE as much as possible.

        Args:
            pubkey (Path): Path to the public key file.
            privkey (Path): Path to the private key file.
            enc (Path): Path to the encrypted data file.
            submission (Path): Path to the submitted file.
        """
        he_context = pyhelayers.SealCkksContext()
        he_context.load_from_file(str(pubkey))
        he_context.load_secret_key_from_file(str(privkey))
        encoder = pyhelayers.Encoder(he_context)

        cx = encoder.encode_encrypt([])
        cx.load_from_file(str(enc))
        x = encoder.decrypt_decode_double(cx)
        cy_sub = encoder.encode_encrypt([])
        cy_sub.load_from_file(str(submission))
        y_sub = encoder.decrypt_decode_double(cy_sub)

        y_true = np.maximum(x, 0.0)
        mae = np.mean(np.abs(y_true - y_sub))  # less is better
        return float(mae)

早速HECCONに取り掛かろうとしたのですが、この問題で使用されているpyhelayersというライブラリをインポートしようとしたところIllegal Instructionと表示されてしまいました。 原因は仮想マシン上で動かそうとしたことなので、ホスト上で直接動かして解決しましたが、このトラブルで割と時間を取られてしまいました。

気を取り直してchallenge.pyを確認したところ、入力 x  -3 \le x \le 3 を満たす一様乱数であったので、とりあえず 0 を暗号化したものを提出しました (Round: 20, MAE: 0.7354, Rank: 1)。

  • heccon_stage_1_calculate_1.py
#!/usr/bin/env python3

# Round: 20, MAE: 0.7354, Rank: 1

import numpy as np
import pyhelayers


N = 2**13

he_context = pyhelayers.SealCkksContext()
he_context.load_from_file('pubkey')
encoder = pyhelayers.Encoder(he_context)

y = encoder.encode_encrypt(np.zeros(N))

y.save_to_file('enc_y')

実はこのときはまだ出力 y = f \! \left( x \right) の計算の存在に気づいていなかったのですが、challenge.pyを見返したところ y = \max \left\{ x, 0 \right\} となっていることに気づきました。 また、 x を暗号化した結果があるので、それを提出すれば何かしらの誤差の評価値が付くことにも気づきました。 逆にいえばここまでそのようなことをしていなかったのでこれまでのスコアが低くなってしまっていました。

暗号化する値を 0 から他の値に変更してもMAEが改善しないので、他の方法を考えました。  y = \max \left\{ x, 0 \right\} を2次関数で近似したら誤差を減らせそうだと考えました。 試行錯誤で2次関数での近似を行ったところ誤差を減らすことができました (Round: 34, MAE: 0.1369, Rank: 3)。 パラメータを調整して y = \frac{3}{4.5^ 2 - 1.5^ 2} \left( \left( x + 1.5 \right)^ 2 - \frac{1.5^ 2}{2} \right) という2次関数を使用したら、さらに誤差が少なくなりました (Round: 35, MAE: 0.1137, Rank: 4)。

  • heccon_stage_1_calculate_2.py
#!/usr/bin/env python3

# Round: 35, MAE: 0.1137, Rank: 4

import numpy as np
import pyhelayers


N = 2**13

he_context = pyhelayers.SealCkksContext()
he_context.load_from_file('pubkey')
encoder = pyhelayers.Encoder(he_context)
x = encoder.encode_encrypt([])
x.load_from_file('enc')

v = pyhelayers.CTile(x)
v.add_scalar(1.5)
v.square()
v.add_scalar(-1.5 ** 2 / 2)
# v.multiply_scalar(3 / (4.5 ** 2 - 1.5 ** 2 / 2))
v.multiply_scalar(3 / (4.5 ** 2 - 1.5 ** 2))
y = v

y.save_to_file('enc_y')

多項式の次数を増やせば誤差を減らせそうだと思いましたが、手作業でパラメータを設定することに限界を感じたのでパラメータを計算する方法を考えました。 numpypolyfitという関数があったのでそれを使用してパラメータを計算しました (Round: 39, MAE: 0.0444, Rank: 4)。 乗算の回数に制限があるので、多項式の項を増やしすぎると逆に誤差が増えたり計算ができなくなったりしました。 いくつか試した結果、多項式の次数が 8 の場合に最も誤差が小さくなりました (Round: 41, MAE: 0.0185, Rank: 2)。

  • heccon_stage_1_calculate_3.py
#!/usr/bin/env python3

# Round: 41, MAE: 0.0185, Rank: 2

import numpy as np
import pyhelayers


x = np.linspace(-3, 3, 100 + 1)
y = np.maximum(x, 0)
# coeff = np.polyfit(x, y, 4)
coeff = np.polyfit(x, y, 8)
# coeff = [-0.0007202732315543946, 4.6793090224511166e-18, 0.015207025357414545, -6.518917576118662e-17, -0.11599788384384442, 2.0608196344928882e-16, 0.48977000799286685, 0.4999999999999998, 0.10149556573347847]


N = 2**13

he_context = pyhelayers.SealCkksContext()
he_context.load_from_file('pubkey')
encoder = pyhelayers.Encoder(he_context)
x = encoder.encode_encrypt([])
x.load_from_file('enc')

v = pyhelayers.CTile(x)
v.sub(v)
v.add_scalar(coeff[0])
for i in coeff[1:]:
    v.multiply(x)
    v.add_scalar(i)
y = v

y.save_to_file('enc_y')

最終ラウンド(Round 42)もMAEは0.0185で順位は一つ下がって3位でした。 ステージ1全体では前半にスコアを稼げなかったことが影響して、スコアの合計が425 ptで順位は5位となりました。

Stage 2 (Score: 535, Rank: 3)

ステージ2の順位表は以下の通りです。

ステージ2のchallenge.pyは以下の通りです。

  • challenge.py
from pathlib import Path

import numpy as np
import pyhelayers


N = 2**12
M = 2**1


class Challenge:
    @classmethod
    def generate(cls, pubkey: Path, dist_dir: Path):
        """Generate the challenge data.

        dist_dir is a distributed directory to contestants.

        Args:
            pubkey (Path): Path to the public key file.
            dist_dir (Path): Path to the directory where the challenge data is saved.
        """
        he_context = pyhelayers.SealCkksContext()
        he_context.load_from_file(str(pubkey))
        encoder = pyhelayers.Encoder(he_context)
        mu = np.random.uniform(-1, 1, (N, 1))
        sigma = np.random.uniform(0.5, 1.0, (N, 1))
        x = np.random.normal(mu, sigma, (N, M))
        x = x.reshape(-1)
        cx = encoder.encode_encrypt(x)
        cx.save_to_file(str(dist_dir / "enc"))
        (dist_dir / "challenge.py").write_text(Path(__file__).read_text())

    @classmethod
    def evaluate(
        cls, pubkey: Path, privkey: Path, enc: Path, submission: Path
    ) -> float:
        """Evaluate the submission.

        submission is a path to a submitted file and the server evaluate it.
        Contestants are supposed to get better MAE as much as possible.

        Args:
            pubkey (Path): Path to the public key file.
            privkey (Path): Path to the private key file.
            enc (Path): Path to the encrypted data file.
            submission (Path): Path to the submitted file.
        """
        he_context = pyhelayers.SealCkksContext()
        he_context.load_from_file(str(pubkey))
        he_context.load_secret_key_from_file(str(privkey))
        encoder = pyhelayers.Encoder(he_context)

        cx = encoder.encode_encrypt([])
        cx.load_from_file(str(enc))
        x = encoder.decrypt_decode_double(cx)
        x = x.reshape(N, M)
        cy_sub = encoder.encode_encrypt([])
        cy_sub.load_from_file(str(submission))
        y_sub = encoder.decrypt_decode_double(cy_sub)
        y_sub = y_sub.reshape(N, M)[:, 0]

        y_true = x.max(axis=1)
        mae = np.mean(np.abs(y_true - y_sub))  # less is better
        return float(mae)

前のステージで「暗号化された入力 x をそのまま提出する」を覚えたので、最初から誤差の評価値が得られて最下位を回避することができました (Round: 43, MAE: 0.4124, Rank: 1)。

このステージで計算するものは2つの数の最大値です。 もう少し詳しくいうと、HECCONの入力 x は全てのステージ共通で 2^ {13} の成分をもつベクトルとして与えられます。(出力 y も同様です。)

 \begin{align}
x &= \begin{bmatrix}
x _0 & x _ 1 & x _ 2 & x _ 3 & \cdots & x _ {2^ {13} - 2} & x _ {2^ {13} - 1}
\end{bmatrix} \\
y &= \begin{bmatrix}
y _0 & y _ 1 & y _ 2 & y _ 3 & \cdots & y _ {2^ {13} - 2} & y _ {2^ {13} - 1}
\end{bmatrix}
\end{align}

そのベクトルの偶数番目の成分と奇数番目の成分の最大値を計算する必要があります。

 \begin{align}
y _ {2 k} &= \max \left\{ x _ {2 k}, x _ {2 k + 1} \right\}
\end{align}

まず最初に y  0 より大きい適当な数で埋めて最大値を雑に近似できないかと考えましたが、実際に試してみても誤差は改善しませんでした。

次に2つの数の平均を計算して、平均より少し大きい数として最大値を近似できないかと考えました。 ただし、この計算をする際に1つ問題があります。 暗号化したままできる演算として加算と乗算がありますがどちらもベクトルの同じ位置の成分に対する演算であり、 0 番目と 1 番目の成分を使って計算するというようなことができないです。

何か方法はないかとHElayersのドキュメントを確認していたのですが、「FHE overview」というページの「What are low level primitives?」という項目に「Low level primitives are the base operations that can be performed directly on ciphertext and include multiplication, addition, and rotation.」と記載されているのを見つけたので、加算と乗算以外にrotationという操作がありそうだということが分かりました。

さらにHElayersのAPIリファレンスを確認したら、「rotate」という関数がありました。 rotate関数は指定した数だけベクトルの成分を前にずらす関数で、例えばベクトル a

 \begin{align}
a &= \begin{bmatrix}
a _0 & a _ 1 & a _ 2 & a _ 3 & \cdots & a _ {2^ {13} - 2} & a _ {2^ {13} - 1}
\end{bmatrix}
\end{align}

としたときに、 \mathrm{rotate} \! \left( a, 2 \right) は以下のようになります。

 \begin{align}
\mathrm{rotate} \! \left( a, 2 \right) &= \begin{bmatrix}
a _ 2 & a _ 3 & \cdots & a _ {2^ {13} - 2} & a _ {2^ {13} - 1} & a _0 & a _ 1
\end{bmatrix}
\end{align}

rotate関数により平均を求められるようになったので、平均より少し大きい数として最大値を近似して誤差を小さくすることができました (Round: 51, MAE: 0.2832, Rank: 3)。 さらに平均に足す値を調整して誤差をさらに改善しました (Round: 52, MAE: 0.27, Rank: 3)。

  • heccon_stage_2_calculate_1.py
#!/usr/bin/env python3

# Round: 52, MAE: 0.27, Rank: 3

import numpy as np
import pyhelayers


N = 2**12
M = 2**1

he_context = pyhelayers.SealCkksContext()
he_context.load_from_file('pubkey')
encoder = pyhelayers.Encoder(he_context)
x = encoder.encode_encrypt([])
x.load_from_file('enc')

v1 = pyhelayers.CTile(x)
v2 = pyhelayers.CTile(x)
v2.rotate(1)
v1.add(v2)
v1.multiply_scalar(1 / 2)
# v1.add_scalar(1 / 2)
v1.add_scalar(1 / 4)
y = v1

y.save_to_file('enc_y')

さらに誤差を減らすために平均に足す値を固定値ではなくデータによって変えることができないかと考えて、 a  b の最大値を求めるために以下の式を使用すれば良いことに気づきました。

 \begin{align}
\max \left\{ a, b \right\} &= \frac{a + b}{2} + \frac{\sqrt{ \left( a - b \right)^ 2 }}{2}
\end{align}

平方根の計算を多項式で近似してこの式を計算することにより誤差を減らすことができました (Round: 57, MAE: 0.0869, Rank: 2)。 関数を多項式で近似する際に変数が取る値の範囲を変えると得られる多項式も変化しますが、その変数の範囲を調整することで誤差をさらに削減できました (Round: 63, MAE: 0.0639, Rank: 2)。

  • heccon_stage_2_calculate_2.py
#!/usr/bin/env python3

# Round: 63, MAE: 0.0639, Rank: 2

import numpy as np
import pyhelayers


# x = np.linspace(0, 4, 100 + 1)
x = np.linspace(0, 6.5, 100 + 1)
y = np.sqrt(x)
coeff = np.polyfit(x, y, 2)
# coeff = [-0.0360415426859049, 0.5508006545032299, 0.4149726744538779]


N = 2**12
M = 2**1

he_context = pyhelayers.SealCkksContext()
he_context.load_from_file('pubkey')
encoder = pyhelayers.Encoder(he_context)
x = encoder.encode_encrypt([])
x.load_from_file('enc')

v1 = pyhelayers.CTile(x)
v2 = pyhelayers.CTile(x)
v3 = pyhelayers.CTile(x)
v4 = pyhelayers.CTile(x)
v3.rotate(1)
v1.add(v3)
v2.sub(v3)
v2.square()

v4.sub(v4)
v4.add_scalar(coeff[0])
for i in coeff[1:]:
    v4.multiply(v2)
    v4.add_scalar(i)

v1.multiply_scalar(1 / 2)
v4.multiply_scalar(1 / 2)
v1.add(v4)
y = v1

y.save_to_file('enc_y')

先ほどは2乗してから平方根を計算してていましたが、絶対値を直接計算すれば良いことに気づきました。

 \begin{align}
\max \left\{ a, b \right\} &= \frac{a + b}{2} + \frac{ \left\lvert a - b \right\rvert }{2}
\end{align}

計算している内容はそれほど変わらないですが、なぜかこの式では多項式の次数をより大きくすることができて誤差を減らすことができました (Round: 69, MAE: 0.05, Rank: 4)。 その後にパラメータを調整してさらに誤差を減らしました (Round: 81, MAE: 0.0361, Rank: 4)。

  • heccon_stage_2_calculate_3.py
#!/usr/bin/env python3

# Round: 81, MAE: 0.0361, Rank: 4

import numpy as np
import pyhelayers


# x = np.linspace(-5, 5, 100 + 1)
x = np.linspace(-3.75, 3.75, 100 + 1)
y = np.abs(x)
coeff = np.polyfit(x, y, 8)
# coeff = [-0.00030210448962013933, -6.133263921907128e-19, 0.0099660761382349, 3.420776859274405e-17, -0.11878183305609469, -2.586859724808885e-16, 0.7836320127885834, 2.546427059557744e-16, 0.2537389143336933]


N = 2**12
M = 2**1

he_context = pyhelayers.SealCkksContext()
he_context.load_from_file('pubkey')
encoder = pyhelayers.Encoder(he_context)
x = encoder.encode_encrypt([])
x.load_from_file('enc')

v1 = pyhelayers.CTile(x)
v2 = pyhelayers.CTile(x)
v3 = pyhelayers.CTile(x)
v4 = pyhelayers.CTile(x)
v3.rotate(1)
v1.add(v3)
v2.sub(v3)

v4.sub(v4)
v4.add_scalar(coeff[0])
for i in coeff[1:]:
    v4.multiply(v2)
    v4.add_scalar(i)

v1.multiply_scalar(1 / 2)
v4.multiply_scalar(1 / 2)
v1.add(v4)
y = v1

y.save_to_file('enc_y')

最終ラウンド(Round 84)も誤差と順位は変わらず、MAEは0.0361で4位でした。 ステージ2全体ではスコアの合計が535 ptで順位は3位となりました。

Stage 3 (Score: 555, Rank: 3)

ステージ3の順位表は以下の通りです。

ステージ3のchallenge.pyは以下の通りです。

  • challenge.py
from pathlib import Path

import numpy as np
import pyhelayers


N = 2**8
M = 2**5
EPS = 1e-5


class Challenge:
    @classmethod
    def generate(cls, pubkey: Path, dist_dir: Path):
        """Generate the challenge data.

        dist_dir is a distributed directory to contestants.

        Args:
            pubkey (Path): Path to the public key file.
            dist_dir (Path): Path to the directory where the challenge data is saved.
        """
        he_context = pyhelayers.SealCkksContext()
        he_context.load_from_file(str(pubkey))
        encoder = pyhelayers.Encoder(he_context)
        mu = np.random.uniform(-1, 1, (N, 1))
        sigma = np.random.uniform(0.5, 1.0, (N, 1))
        x = np.random.normal(mu, sigma, (N, M))
        x = x.reshape(-1)
        cx = encoder.encode_encrypt(x)
        cx.save_to_file(str(dist_dir / "enc"))
        (dist_dir / "challenge.py").write_text(Path(__file__).read_text())

    @classmethod
    def evaluate(
        cls, pubkey: Path, privkey: Path, enc: Path, submission: Path
    ) -> float:
        """Evaluate the submission.

        submission is a path to a submitted file and the server evaluate it.
        Contestants are supposed to get better MAE as much as possible.

        Args:
            pubkey (Path): Path to the public key file.
            privkey (Path): Path to the private key file.
            enc (Path): Path to the encrypted data file.
            submission (Path): Path to the submitted file.
        """
        he_context = pyhelayers.SealCkksContext()
        he_context.load_from_file(str(pubkey))
        he_context.load_secret_key_from_file(str(privkey))
        encoder = pyhelayers.Encoder(he_context)

        cx = encoder.encode_encrypt([])
        cx.load_from_file(str(enc))
        x = encoder.decrypt_decode_double(cx)
        x = x.reshape(N, M)
        cy_sub = encoder.encode_encrypt([])
        cy_sub.load_from_file(str(submission))
        y_sub = encoder.decrypt_decode_double(cy_sub)
        y_sub = y_sub.reshape(N, M)

        mu = x.mean(axis=1, keepdims=True)
        sigma2 = x.var(axis=1, keepdims=True)
        y_true = (x - mu) / np.sqrt(sigma2 + EPS)
        mae = np.mean(np.abs(y_true - y_sub))  # less is better
        return float(mae)

このステージでも最初のラウンドでは暗号化された入力 x をそのまま提出しました (Round: 85, MAE: 0.549, Rank: 1)。

このステージで行う計算は y = \frac{x - \mu}{\sqrt{\sigma^ 2 + \epsilon}} というものです。 ただし、 \mu, \sigma^ 2 はそれぞれ M = 2^ 5 個ごとに分けたデータの平均と分散であり、 \epsilon = 10^ {-5} です。 ベクトルのどの成分に対して行う計算なのかを明示すると以下のような式になります。

 \begin{align}
\mu _ {i M} &= \frac{1}{M} \sum _ {j = 0}^ {M - 1} x _ {i M + j} \\
\sigma^ 2 _ {i M} &= \frac{1}{M} \sum _ {j = 0}^ {M - 1} \left( x _ {i M + j} - \mu _ {i M} \right)^ 2 \\
y _ {i M + j} &= \frac{x _ {i M + j} - \mu _ {i M}}{\sqrt{\sigma^ 2 _ {i M} + \epsilon}}
\end{align}

平均 \mu の計算は基本的にはステージ2と同じですが、普通に計算すると m = 5 としたときに M = 2^ m 回の加算が必要になります。 そこで、ベクトル a に対して、

 \begin{align}
b &= a + \mathrm{rotate} \! \left( a, 1 \right) \\
c &= b + \mathrm{rotate} \! \left( b, 2 \right) \\
d &= c + \mathrm{rotate} \! \left( c, 4 \right) \\
&\;\;\vdots
\end{align}

という計算をすると、それぞれのベクトルの成分は

 \begin{align}
b _ i &= a _ i + a _ {i + 1} \\
c _ i &= b _ i + b _ {i + 2} = a _ i + a _ {i + 1} + a _ {i + 2} + a _ {i + 3} \\
d _ i &= c _ i + c _ {i + 4} = a _ i + a _ {i + 1} + a _ {i + 2} + a _ {i + 3} + a _ {i + 4} + a _ {i + 5} + a _ {i + 6} + a _ {i + 7} \\
&\;\;\vdots
\end{align}

となるので、 m 回の加算で M = 2^ m 個の値の和を求めることができます。

加算の回数は恐らく誤差には影響しないですが、計算時間が短縮されるので試行錯誤の回数を増やせました。 また、手元の環境で平均絶対誤差の評価を行う方法を1日目終了後の夜間の時間に確認したので、ステージ3からは試行錯誤の回数をさらに増やすことができました*2

分散 \sigma^ 2 の計算は、以下の式を利用して計算しました。

 \begin{align}
\sigma^ 2 &= \frac{1}{M} \sum _ {j = 0}^ {M - 1} x _ j^ 2 - \left( \frac{1}{M} \sum _ {j = 0}^ {M - 1} x _ j \right)^ 2 \\
&= \frac{1}{M} \sum _ {j = 0}^ {M - 1} x _ j^ 2 - \mu^ 2
\end{align}

 \sigma^ 2 が計算できたら、 f \! \left( t \right) = \frac{1}{\sqrt{t + \epsilon}} 多項式近似して y = \left( x - \mu \right) f \! \left( \sigma^ 2 \right) というようにして y を計算します。 ただし、 \mu, \sigma^ 2 の計算結果が格納されているベクトルは、 i M 番目の成分に必要な値があり他の成分には不要な値が入ってしまっています。 そこで、以下のベクトル m (プログラム中ではmask)を用意して乗算することにより必要な部分のみを取り出して計算するようにしました。

 \begin{align}
m &= \left[ \begin{array}{cccc|cccc|c|cccc}
1 & 0 & \cdots & 0 & 1 & 0 & \cdots & 0 & \cdots & 1 & 0 & \cdots & 0
\end{array} \right]
\end{align}

このようにして y を計算することにより誤差を減らすことができました (Round: 108, MAE: 0.0168, Rank: 1)。 そしてパラメータを調整してさらに誤差を減らしました (Round: 112, MAE: 0.0012, Rank: 1)。 この辺りのラウンドは、私たちのチームが1位でそのすぐ後ろに2位のチームが追いかけているような展開だったのでとても緊張した記憶があります。

  • heccon_stage_3_calculate_1.py
#!/usr/bin/env python3

# Round: 112, MAE: 0.0012, Rank: 1

import numpy as np
import pyhelayers


EPS = 1e-5

# x = np.linspace(0.5, 3.25, 100 + 1)
x = np.linspace(0.25, 1.5, 100 + 1)
y = 1 / np.sqrt(x + EPS)
coeff = np.polyfit(x, y, 6)
# coeff = [1.5857495139150735, -9.607447497530757, 24.01622359716313, -32.07290856278935, 24.834515682645815, -11.489790140032953, 3.732927571909671]


n, m = 8, 5
N = 2 ** n
M = 2 ** m

he_context = pyhelayers.SealCkksContext()
he_context.load_from_file('pubkey')
encoder = pyhelayers.Encoder(he_context)
x = encoder.encode_encrypt([])
x.load_from_file('enc')

v2 = pyhelayers.CTile(x)
for i in range(m):
    v3 = pyhelayers.CTile(v2)
    v3.rotate(1 << i)
    v2.add(v3)
v2.multiply_scalar(1 / M)
mu_ = v2

v6 = pyhelayers.CTile(x)
v6.square()
for i in range(m):
    v7 = pyhelayers.CTile(v6)
    v7.rotate(1 << i)
    v6.add(v7)
v6.multiply_scalar(1 / M)
v11 = pyhelayers.CTile(mu_)
v11.square()
v6.sub(v11)
var_ = v6

v8 = encoder.encode_encrypt(np.zeros(N * M))
v8.add_scalar(coeff[0])
for i in coeff[1:]:
    v8.multiply(var_)
    v8.add_scalar(i)
inv_sigma_ = v8

v4 = np.zeros((N, M))
v4[:, 0] = 1
v4 = encoder.encode_encrypt(v4.reshape(-1))
mask = v4

v10 = encoder.encode_encrypt(np.zeros(N * M))
for i in range(M):
    v12 = pyhelayers.CTile(x)
    v12.rotate(i)
    v12.sub(mu_)
    v12.multiply(inv_sigma_)
    v12.multiply(mask)
    v12.rotate(-i)
    v10.add(v12)
y = v10

y.save_to_file('enc_y')

numpypolyfit関数は最小二乗法で近似多項式を求めていますが、この問題では平均絶対誤差(MAE)によって評価されるので最小絶対偏差法で近似多項式を求めた方がMAEが小さくなるのではないかと考えました。 最小絶対偏差法の計算はIteratively Reweighted Least Squares(IRLS)というアルゴリズムを使用して行いました。 最小絶対偏差法で求めた近似多項式を使用したところ、さらに誤差が改善されて1位を奪還することができました (Round: 118, MAE: 0.0009, Rank: 1)*3

ただし、CTF終了後に別のデータで試したときは誤差が改善しなかったので、最小絶対偏差法を使用したことによる効果で誤差が改善したのかどうかは微妙です。 効果があるのかは不明ですが、この問題はある程度方針が固まったらあとはパラメータを調整して最適化する作業になるので、その際に試せる方法は多い方が良いとは思います。

  • heccon_stage_3_calculate_2.py
#!/usr/bin/env python3

# Round: 118, MAE: 0.0009, Rank: 1

import numpy as np
import pyhelayers


EPS = 1e-5

x = np.linspace(0.2, 1.3, 1000 + 1)
y = 1 / np.sqrt(x + EPS)

X = np.array([x ** i for i in range(6 + 1)]).T
coeff = np.zeros(6 + 1)
e = y - X @ coeff
for _ in range(1000):
    w = 1 / np.sqrt(np.maximum(np.abs(e), np.finfo(float).eps))
    coeff, e, *_ = np.linalg.lstsq(w[:, np.newaxis] * X, w * y)
coeff = list(reversed(coeff))
# coeff = [4.56120920586022, -23.65861092777583, 50.54141914087366, -57.55204781957862, 37.88955927295624, -14.851971699268269, 4.071053412365981]


n, m = 8, 5
N = 2 ** n
M = 2 ** m

he_context = pyhelayers.SealCkksContext()
he_context.load_from_file('pubkey')
encoder = pyhelayers.Encoder(he_context)
x = encoder.encode_encrypt([])
x.load_from_file('enc')

v2 = pyhelayers.CTile(x)
for i in range(m):
    v3 = pyhelayers.CTile(v2)
    v3.rotate(1 << i)
    v2.add(v3)
v2.multiply_scalar(1 / M)
mu_ = v2

v6 = pyhelayers.CTile(x)
v6.square()
for i in range(m):
    v7 = pyhelayers.CTile(v6)
    v7.rotate(1 << i)
    v6.add(v7)
v6.multiply_scalar(1 / M)
v11 = pyhelayers.CTile(mu_)
v11.square()
v6.sub(v11)
var_ = v6

v8 = encoder.encode_encrypt(np.zeros(N * M))
v8.add_scalar(coeff[0])
for i in coeff[1:]:
    v8.multiply(var_)
    v8.add_scalar(i)
inv_sigma_ = v8

v4 = np.zeros((N, M))
v4[:, 0] = 1
v4 = encoder.encode_encrypt(v4.reshape(-1))
mask = v4

v10 = encoder.encode_encrypt(np.zeros(N * M))
for i in range(M):
    v12 = pyhelayers.CTile(x)
    v12.rotate(i)
    v12.sub(mu_)
    v12.multiply(inv_sigma_)
    v12.multiply(mask)
    v12.rotate(-i)
    v10.add(v12)
y = v10

y.save_to_file('enc_y')

この後も誤差を減らせないかと試行錯誤しましたが0.0009が最も良い値で、最終ラウンド(Round 126)までそのまま逃げ切り最終ラウンドの順位は1位でした。 ステージ3全体では前半にプログラムの実装に苦戦したことが影響して、スコアの合計が555 ptで順位は3位となりました。

Stage 4 (Score: 607, Rank: 2)

ステージ4の順位表は以下の通りです。

ステージ4のchallenge.pyは以下の通りです。

  • challenge.py
from pathlib import Path

import numpy as np
import pyhelayers


N = 2**9
M = 2**4


class Challenge:
    @classmethod
    def generate(cls, pubkey: Path, dist_dir: Path):
        """Generate the challenge data.

        dist_dir is a distributed directory to contestants.

        Args:
            pubkey (Path): Path to the public key file.
            dist_dir (Path): Path to the directory where the challenge data is saved.
        """
        he_context = pyhelayers.SealCkksContext()
        he_context.load_from_file(str(pubkey))
        encoder = pyhelayers.Encoder(he_context)
        mu = np.random.uniform(-1, 1, (N, 1))
        sigma = np.random.uniform(0.5, 1.0, (N, 1))
        x = np.random.normal(mu, sigma, (N, M))
        x = x.reshape(-1)
        cx = encoder.encode_encrypt(x)
        cx.save_to_file(str(dist_dir / "enc"))
        (dist_dir / "challenge.py").write_text(Path(__file__).read_text())

    @classmethod
    def evaluate(
        cls, pubkey: Path, privkey: Path, enc: Path, submission: Path
    ) -> float:
        """Evaluate the submission.

        submission is a path to a submitted file and the server evaluate it.
        Contestants are supposed to get better MAE as much as possible.

        Args:
            pubkey (Path): Path to the public key file.
            privkey (Path): Path to the private key file.
            enc (Path): Path to the encrypted data file.
            submission (Path): Path to the submitted file.
        """
        he_context = pyhelayers.SealCkksContext()
        he_context.load_from_file(str(pubkey))
        he_context.load_secret_key_from_file(str(privkey))
        encoder = pyhelayers.Encoder(he_context)

        cx = encoder.encode_encrypt([])
        cx.load_from_file(str(enc))
        x = encoder.decrypt_decode_double(cx)
        x = x.reshape(N, M)
        cy_sub = encoder.encode_encrypt([])
        cy_sub.load_from_file(str(submission))
        y_sub = encoder.decrypt_decode_double(cy_sub)
        y_sub = y_sub.reshape(N, M)[:, 0]

        y_true = x.max(axis=1)
        mae = np.mean(np.abs(y_true - y_sub))  # less is better
        return float(mae)

このステージでも最初のラウンドでは暗号化された入力 x をそのまま提出しました (Round: 127, MAE: 1.3694, Rank: 1)。

このステージでは、 M = 2^ 4 個の値の最大値を計算します。

 \begin{align}
y _ {i M} &= \max \left\{ x _ {i M}, x _ {i M + 1}, \dots, x _ {i M + M - 1} \right\}
\end{align}

これまでのステージでは前半の順位が低くステージ全体のスコアが伸びなかった面があったので、このステージでは最初から誤差を減らして高い順位を取りたいと考えました。 ステージ4でも、ステージ2のように平均より少し大きい数として最大値を近似できないかと考えて試してみました。 平均に足した 1.5 という値が丁度良かったのか、ステージ4の2回目のラウンドも1位になりました (Round: 128, MAE: 0.3861, Rank: 1)。 この後さらにパラメータを調整して誤差を改善しました (Round: 130, MAE: 0.3498, Rank: 2)。

  • heccon_stage_4_calculate_1.py
#!/usr/bin/env python3

# Round: 130, MAE: 0.3498, Rank: 2

import numpy as np
import pyhelayers


n, m = 9, 4
N = 2 ** n
M = 2 ** m

he_context = pyhelayers.SealCkksContext()
he_context.load_from_file('pubkey')
encoder = pyhelayers.Encoder(he_context)
x = encoder.encode_encrypt([])
x.load_from_file('enc')

v2 = pyhelayers.CTile(x)
for i in range(m):
    v3 = pyhelayers.CTile(v2)
    v3.rotate(1 << i)
    v2.add(v3)
v2.multiply_scalar(1 / M)
mu = v2

v3 = pyhelayers.CTile(mu)
# v3.add_scalar(1.5)
v3.add_scalar(1.3)
y = v3

y.save_to_file('enc_y')

さらに誤差を減らすための方法を考えて、標準偏差 \sigma を計算すれば最大値と平均の差の目安になるのではないかということを思いつきました。 具体的には平均 \mu 標準偏差 \sigma を使用して \mu + \alpha \sigma というようにして最大値を近似できるのではないかと考えました。 分散 \sigma^ 2 の計算はステージ3で行ったので、平方根多項式で近似すれば標準偏差 \sigma の計算ができます。 あとは最大値と平均の差が標準偏差 \sigma の何倍かという値 \alpha も必要ですが、データの個数 M = 2^ 4 = 16 を見てとりあえず決め打ちで \alpha = 1.8 としました。

この 1.8 という値がほぼ最適であったのか、この計算結果を提出したラウンドでも1位になりました (Round: 135, MAE: 0.2197, Rank: 1)。 この計算についてもパラメータを調整して誤差をさらに改善しました (Round: 150, MAE: 0.2191, Rank: 2)。

  • heccon_stage_4_calculate_2.py
#!/usr/bin/env python3

# Round: 148, MAE: 0.2192, Rank: 2

import numpy as np
import pyhelayers


x = np.linspace(0, 2.5, 1000 + 1)
y = np.sqrt(x)

X = np.array([x ** i for i in range(7 + 1)]).T
coeff = np.zeros(7 + 1)
e = y - X @ coeff
for _ in range(1000):
    w = 1 / np.sqrt(np.maximum(np.abs(e), np.finfo(float).eps))
    coeff, e, *_ = np.linalg.lstsq(w[:, np.newaxis] * X, w * y)
coeff = list(reversed(coeff))
# coeff = [0.08507371914473838, -0.8199410426728595, 3.214640799893507, -6.603337834146714, 7.666560640974134, -5.183044184431915, 2.546094604486351, 0.09553325793210383]


n, m = 9, 4
N = 2 ** n
M = 2 ** m

he_context = pyhelayers.SealCkksContext()
he_context.load_from_file('pubkey')
encoder = pyhelayers.Encoder(he_context)
x = encoder.encode_encrypt([])
x.load_from_file('enc')

v2 = pyhelayers.CTile(x)
for i in range(m):
    v3 = pyhelayers.CTile(v2)
    v3.rotate(1 << i)
    v2.add(v3)
v2.multiply_scalar(1 / M)
mu = v2

v6 = pyhelayers.CTile(x)
v6.square()
for i in range(m):
    v4 = pyhelayers.CTile(v6)
    v4.rotate(1 << i)
    v6.add(v4)
v6.multiply_scalar(1 / M)
v5 = pyhelayers.CTile(mu)
v5.square()
v6.sub(v5)
var = v6

v8 = encoder.encode_encrypt(np.zeros(N * M))
v8.add_scalar(coeff[0])
for i in coeff[1:]:
    v8.multiply(var)
    v8.add_scalar(i)
sigma = v8

v3 = pyhelayers.CTile(mu)
v7 = pyhelayers.CTile(sigma)
v7.multiply_scalar(1.8)
v3.add(v7)
y = v3

y.save_to_file('enc_y')

この後も別の方法を試したりしましたが、これ以上誤差を減らすことはできなかったです*4。 Writeup執筆時点では他にも方法があったのではないかと思うのですが、CTF中はこれ以上の方法は思いつかなかったです。

最終ラウンド(Round 168)までMAEは0.2191のままで、最終ラウンドの順位は4位でした。 ステージ4全体では前半に良い順位が取れたので、スコアの合計が607 ptで順位は2位となりました。

*1: Writeupを書くまでこのファイルの存在に気づいておらず、HECCON開始時はHElayersのAPIリファレンスを睨みつけながらプログラムを書いていました。

*2: 逆に言えば1日目はそこまでする余裕がなかったです。

*3: MAEをもう1桁小さくすれば表示上MAEが0になって1位を取り続けることができますね

*4: 見当違いなことをしてMAEが20000ぐらいになってしまったデータもありました。