DEIM Forum 2023 4a-6-3

周期解析による変動天体の検出
千原 直己†

高田唯史††

藤原 靖宏†††

鬼塚真†† 大阪大学大学院情報科学研究科 〒 565-0871 大阪府吹田市山田丘 1–5††† NTT コミュニケーション科学基礎研究所 〒 243-0198 神奈川県厚木市森の里若宮 3-1†† 自然科学研究機構国立天文台 〒 181-8588 東京都三鷹市大沢 2-21-1†† 総合研究大学院大学 〒 181-8588 東京都三鷹市大沢 2-21-1
E-mail: †{chihara.naoki,onizuka}@ist.osaka-u.ac.jp, ††tadafumi.takata@nao.ac.jp,†††yasuhiro.fujiwara.kh@hco.ntt.co.jp
あらまし

変動天体とは, 光度や位置などが時間とともに変化する天体のことであり, その検出は, 高精度な情報や新たな知見の獲得に繋がるため重要である. 既存技術では, 天体ごとの統計量ベースの探索手法が主流であるが, 問題が2 つ存在する. i) 観測情報を統計的な数値に集約してしまうため, 変動天体の持つとされている周期性関連の情報が損なわれてしまう, ii) 天体の観測可能な時刻が限定されることにより, 基本的に多数の欠損値が含まれたデータしか用意できないことである. これらの問題を解決するため, 本稿では, スパースモデリングを取り入れた周期解析による変動
天体の検出を提案する. 具体的には, 周期解析により時系列の重要な情報を抽出した後, それらを活用し, 二値分類にて変動天体かどうかの判定を行う. 実験では, ノイズの混入した人工データに対しても正常な学習が可能であることを明確にした. 加えて, 実データに対しても同様に提案手法の有効性を示した.
キーワード 変動天体, 周期解析, スパースモデリング, 異常検知
損値への耐性が低い, ii) 時系列に複数の周期的な変動が含まれ1 はじめに
変動天体には, i) 位置が変動する天体, ii) 光度が変動する天体の 2 種類存在する. 変動については, 長期的なものもあれば,局所的なもの, それらが混在しているものに加え, 一過性であるものや周期性を持つものもなど多岐にわたる [1]. 加えて, それは全天体の約 1% 程度を占める. 変動天体の代表例としては,
閃光星, 新星, 超新星などが挙げられる. 本研究では, 光度が周期的に変動する変動天体のみを対象とする. 非常に稀な天体現
象を持つ変動天体を発見することは, 宇宙の構造, 進化に関する研究や, 物理現象の解明などにとって有益となりうる高精度な情報や新たな知見を得ることと密接な関わりがある [2–4].
既存手法には統計的手法を活用したものが多い. ここでは有名な手法を二つ紹介する. 一つ目は B.Sesar らの手法 [2] である. この手法は, 時系列から得られる intrinsic variability σ,

χ2 , 等級の測定値の平均, 標準偏差, 歪度といった統計的数値を活用したスコアリングを基にした変動天体検出法の一種である. この手法では, 天体データを統計的な数値に集約するため, 変動天体の持つ時間情報が損なわれる, すなわち, 変動天体の特性として知られる光度曲線の周期性を活用出来ていな
い, という問題点が考えられる. 二つ目の手法は, lomb-scargle

periodogram [5, 6] である. こちらは, 最小二乗フィッティングを活用した, 時系列に存在する周期性の分析を行うための手法である. サンプリング間隔が一定でない時系列に対しても適用できるため, 観測環境により一定間隔でのデータ収集が困難である天文分野にてよく活用されている手法である. しかし, この手法にもいくつか問題点がある. 具体的には, i) ノイズや欠る場合にはそれらを同時に分析することができない場合がある,

iii) 周波数分解能が低い, といったものが挙げられる.
本研究では, 機械学習的アプローチを活用した周期解析により既存手法の問題点を解決することで, より高精度な変動天体検出の達成を目的とする. 全体の流れとしては, 初めに回帰分析を活用することで天体データから周期情報を取得する. そして, 後半部分では, 得た情報を元に特徴量設計を行い, 分類タスクの解決により変動天体候補を取得する. しかし, この実現のために幾つかの課題点が存在する. 一つ目が, 利用する天体デー
タに多数の欠損値が含まれていることである (図 1 参照). 原因として考えられることは, 観測の際, 天候や季節などといった外的影響に加え, 望遠鏡の口径による観測範囲や天体の観測時間に関する制限などが挙げられる. そのため, ナイーブな回帰分析を適用すると過学習が生じてしまい, 望ましい結果が得られる見込みが薄い. 二つ目が, 特徴量設計に使用するデータ候
補が周期情報のみであるということである. すなわち, 分類タ
スクによる変動天体の予測が周期解析結果に完全に依存する.
故に, 懸念事項として, 誤った周期解析結果が得られた場合, i)
対象の天体データが含まれたデータセットで訓練を行うことにより, モデルのパフォーマンスが大幅に下がる, ii) 対象の天体
データを正しく予測することが大凡不可能である, といったことが挙げられる.
本研究では, 上記の前者の問題点を解決すべく, スパースモデ
リングの活用を提案する. スパースモデリングとは, データから有為な性質のみを学習し, 全体から重要な変数のみに限定させる, という変数選択手法の一種である. 訓練データから核と等級

表 2: 記号とその定義
記号

定義

周期解析

時間 [⽇]

図 1: 使用データセット例

n

時系列の長さ各天体が観測された明るさの時間変化を示したデータである. 横軸が時m

周期解析に用いる周波数の候補数間であり, 縦軸が等級を示している. 但し, 縦軸の等級の値は実際の値

N

天体の件数

を標準化したものである. 観測されたデータの割合が全体の約 0.01%

fn

ナイキスト周波数と非常にスパースなデータとなっている.

X

説明変数 (周波数) X ∈ Rn×2my

目的変数 (天体データ) y ∈ Rn

Y

天体データセット Y = {y1 , ..., yN }w

回帰係数 w ∈ R2m

ν

周波数候補 ν ∈ Rm

なる箇所を選定するため, ノイズなどの影響による細かなデー
タに対して過学習しにくいといったメリットが挙げられる. また, 後者の問題点を解決するために, ドメイン知識を駆使した特Group Lasso

徴量エンジニアリングを提案する. 分類タスクの訓練及び予測

J

グループ数

を行う際に, 変動天体の傾向を正確に反映させることで高精度

gj

j 番目のグループに属する変数の数

(1 <
=j<
= J)

な検出を実現する.

G

グループセット G = {G1 , ..., GJ }(Gj ∈ Rgj )二値分類

以下に提案手法の貢献点を要約する.

F

特徴量リスト F ∈ RN ×(m+5)

高精度 既存の変動天体検出手法とは異なり, 提案手法では重Fd

ドメイン特徴量リスト Fd ∈ RN ×4

要周期の抽出により, 変動天体の特徴である周期性を考慮

L

変動天体候補セット

した検出を行う. さらに, 重要周期の抽出は数百件程度の候補から自動探索を行うため, 周期長が不明な変動天体に

但し, F は本研究における提案モデルを表す関数である.

対しても高精度な処理が担保される.

2. 2 技術的課題
ロバスト性 本研究で取り扱うデータセットは, 欠損値が非常に多いことに加え, ノイズが混入している蓋然性が高い. このような特徴を持つデータに対しても, スパースモデリン
グを活用することにより, ノイズや過学習に対してロバス
トな処理を実現した.

本研究の課題は主に 2 種類存在する. 1 つ目が, 解析すべき
データセットに多数の欠損値が含まれている点である (図 1 参照). そのため, 天体データに直接フーリエ変換を利用するためには欠損値補完を適用する必要がある. しかし, 観測データ点の割合は平均 0.01%程度であり, 残りの約 99.99%に欠損値補完

高スケール性 既存の異常検知手法であると, 全体の傾向に基を適用すると, 正確な周期解析が困難であるため, フーリエ変換

づいて例外的な挙動の発見をする必要があるため, 並列化を直接適用することが困難である. これを解決するため, 回帰が難しい傾向にある. しかし, 提案手法では, 天体毎に独立問題として本研究を扱う. 但し, その際に過学習にも気を配るに処理が可能であるため, 100 万天体規模に対しても高ス必要がある. こちらの問題を緩和させるため, 特徴量選択によ

ケールな処理が期待できる.

りモデルの汎化性を上げる技術としてスパースモデリングを用いる.

2 つ目が, 変動天体検出の際に使用する特徴量設計である. 現

2 研究背景状, 特徴量設計に使用するデータ候補が周期情報のみであるた

本章では, 解決すべき問題の定義並びに本論文で必要な前提

め, 高精度な分類を実現するためには, 正確な周期解析が求めら知識の説明を行う. 表 2 に提案時に使用する記号とその定義をれる. しかし, 使用する天体データには多数の欠損値が含まれ記載する.

ているため, 精度には限度があると考えられる. また, 例え周期解析により天体の持つ適切な周期情報が得られたとしても, 分

2. 1 問 題 定 義

類の際にそれらが上手く機能する変動天体の傾向に関する情報

本研究の目的は, 多数の天体データから変動天体を検出することである. 入力は天体の光度の変化を表す時系列データセッ
ト Y = {y1 , ..., yN }, 周波数候補の数 m, ナイキスト周波数

fn , 特徴量リスト F に属する各特徴量の内, ドメイン知識に関がなければ高精度な結果は期待できないため, ドメイン知識を考慮した特徴量設計を行う必要がある.

2. 3 前 提 知 識

連したものに限定したもの Fd であり, 出力は変動天体候補セッ

フーリエ変換: 時系列を時間領域から周波数領域へ変換する

ト L である. 本研究にて取り組む問題を数式化すると (1) 式の手法. 本研究では, 天体データの周期性を調べるために用いる.

ように与えられる.

具体的には, Y の一要素であり, 光度の変化を表現している時
系列 y を正弦波, 余弦波の和で書き表すために使用する. 数式L = F (Y, m, fn , Fd )で書き表すと, 以下の (2) 式のようになる. 但し, t は時刻を表(1)しており, νj は周波数候補の j 番目を表している.

不可分な関係にある. Lasso ではこれらを独立に扱った解析しか行えないが, Group Lasso ではグループ化することにより上

y=

m
∑

記の問題を解消した解析が可能である.

sin(2πtνj ) + cos(2πtνj )

(2)

j=1

LightGBM [9]: 決定木アルゴリズムに基づいた勾配ブース

Lasso [7]: スパースモデリングの一種. L1 正則化項を活用

ティングの教師あり機械学習手法 [10] の一種である. すなわすることにより, 更新を繰り返したとしても不必要な特徴量にち, 比較的浅い決定木を複数用意し, それらを束ねて学習を行対応した回帰係数は 0 のまま維持される, すなわち, 特徴量選

う. 目的変数に応じて、説明変数を分類することが可能である.

択が行われる. 詳細については以下の補題 2 にて示し, 本稿に関連した手法に XGBoost [11] があるが, 最終的には 100 万天

おける Lasso の損失関数は以下の (3) 式のように書き表す.

体規模の分類を行う必要があるため, より計算コストが小さい

補題 1. Lasso を用いて回帰分析を行う際, 説明変数 X ∈ Rn×m ,目的変数 y ∈ Rn , 回帰係数 w−k ∈ Rm−1 に対して, wk が

−nλ < wk < nλ を満たす場合, 値の更新が行われない. 但し
w−k とは, 列ベクトル w から k 番目の要素, すなわち wk を除いたベクトルを表す.

LightGBM を採用した.

3 提案手法
本章では, 変動天体を検出するまでの提案モデルについて説明する. 提案モデルは以下の 2 部から構成されている. 概要に関しては図 3 にてまとめている.

ELasso =証明.

1||y − Xw||22 + λ||w||1
2n

(3)

周期解析
複数の周波数候補を表現する, 正弦波及び余弦波にて構成された説明変数 X 及び, 天体データである目的変数 Y を

付録 1 参照

周期解析の際, 重要度の高い特徴量のみを検出できれば良いため, 特徴量選択が可能な Lasso は本研究にて適切な手段であると考えられる. なぜなら, 特徴量選択は分析の際, 過学習の発生を抑制する効果を持つためである. この仕組みについては, 特用意する. その後, 各天体データに対して Group Lasso を適用し, 回帰係数 w, 決定係数 reg score を算出する. (図 3左側の周期解析部分)
変動天体検出

徴量選択によりデータの根幹を表現するために不要と判断され回帰係数 w を始めとした周期情報に加え, ドメイン知識にた特徴量を 0 とすることによるものである. すなわち, データ

よって構成された特徴量セットを用いて, LightGBM によの細部までの表現を行わないことにより, ノイズなどの影響にり変動天体かどうかの判定を行う. (図 3 右側の変動天体検よる些細な誤差への過剰なフィッティングが抑制される.

出部分)

Group Lasso [8]: Lasso の応用手法. Lasso と異なる特徴

3. 1 周 期 解 析として, 特徴量選択の際, 一つ一つに対してではなく, グルー

欠損値が多数混在したデータセットに対して, 周期情報の取プ単位で行う. すなわち, 事前に関係性のある特徴量群が考え得にフーリエ変換を直接適用するためには欠損値補完が必要不うる場合, それらをグループ化し回帰分析に反映させることで,可欠であるが, 欠損値が多すぎるために, 有用な補完値を得る

Lasso 以上の高精度なモデリングが期待できる. 詳細についてことが非常に困難である. 故に, 回帰問題として周期解析を取は以下の補題 1 にて示し, 本稿における Group Lasso の損失関り扱う. このような周期解析を実行するためには, 初めに説明

数は以下の (4) 式のように書き表す. また, 式中の wGj は j 番

変数を用意する必要がある. 説明変数に異なる周波数である正目のグループ Gj に属する特徴量を表すベクトルである.

弦波並びに余弦波をそれぞれ m ずつ, 合計 2m 個用意する. なお, こちらは加藤らの手法 [12] に基づいている. 具体的には説補題 2. Group Lasso を用いて回帰分析を行う際, 説明変数

明変数は以下のように記述することができる.

X ∈ Rn×m , 目的変数 y ∈ Rn , 回帰係数 w ∈ Rm , グループdef

セット G ∈ RJ に対して, ŵGj (⇐⇒ wGj − η∇f (wGj )) が||ŵGj || > ηλ を満たす場合, 値の更新が行われない.

xij =


 cos(2πti νj )

(1 <
= m)
=j<

 sin(2πt ν )(m + 1 <
=j<
= 2m)

i j

(5)

但し, t, ν は以下のように定義する. すなわち, 表現できる最大

EGroupLasso =∑1||y − Xw||22 + λ||wGj ||
2n
G ∈G

(4)

j

周波数が fn である.

fn
(6)
m
ま た Group Lasso を 適 用 す る た め に, グ ル ー プ セ ッ ト
ti = i,

証明.

付録 2 参照

本研究では, 任意の周波数候補に対して, 正弦波及び余弦波の計 2 つの特徴量を与える. 但し, これらは周期解析である以上,

G を 定 義 す る 必 要 が あ る.

νj = j

最終的に抽出したい情報は重要 周 期 長 で あ る た め , 任 意 の ti , あ る 周 波 数 ν j に 対 し て,

⼊⼒

周期解析

天体データ

Group Lasso

変動天体検出
回帰曲線

出⼒
訓練データ
学習

…

…

特徴量設計

モデル

予測

変動天体候補テストデータ

ドメイン知識図 3: 提案手法全体像

3. 2 変動天体検出

Algorithm 1 PeriodicAnalysis
Input: i) 説明変数: X

14

変動天体検出を実現するためには, 以後教師ラベルとして取ii) 目的変数: y

り扱う変動天体識別子 objclass が, 0(非変動天体) か 1(変動天iii) グループセット: G体) かを識別する分類タスクとして取り扱う必要がある.

Output: i) 回帰係数: w3. 2. 1 特徴量について

ii) 決定係数: reg score

こちらでは二値分類タスクにて実際に使用する特徴量リスト

1: // Group Lasso
2: w = 0 // 初期化

F に関して説明する. 特徴量リスト F はドメイン特徴量リス

3: while w is not converged doトである Fd を構成する 4 種類の特徴量, 回帰係数 w, 決定係

4:ŵ ← w − η∇f (w) // ∇f (w) = X ⊤ (Xw − y)/n数 reg score により構成されている. また, Fd とは, ドメイン知5:

for Gj to G do

識を駆使し変動天体検出を実現するために採択に至った特徴量

6:7:

wGj ← max(0, 1 − ηλ/||ŵGj ||)ŵGj

セットであり, 具体的には平均 m ap40, z apertureﬂux 40 mag,end fordiﬀ ap40 から成る. 以下では各要素について説明する [13].

8: end while
9: reg score = ||Xw − y||2 /||y − y||2w 周期解析により算出した天体データの各周波数候補に対応

10: return w, reg score

した周期関数の係数のベクトル. 各係数の大きさは対応す
る周期長の重要度を表す. 各要素は, 特徴量選択により採

cos(2πti νj ), sin(2πti νj ) の組み合わせで周期を表現する. 言い換えると, cos(2πti νj ) のみが特徴量として採択され, sin(2πti νj )は棄却されるといった結果として望ましくないということである. 以上の条件を満たす結果を得るために, グループセット

G = {{1, m + 1}, {2, m + 2}, ..., {m, 2m}} のように表現する 1 .
3. 1. 1 アルゴリズム
周期解析の擬似コードを Algorithm1 に記載する. 2 行目
で は, 零 ベ ク ト ル に て w の 初 期 化 を 行 う.

択された重要周期のみからなるスパースなベクトルである.

reg score 周期解析により求めた回帰曲線が目的変数をどの程度表現出来ているか表す決定係数. 周期解析の精度は特徴量 w の確からしさに相当するため, 変動天体検出に有効であると考えた.
平均 m ap40

観測点に対応した天体の光度である m ap40 を

3-8 行 目 で は,平均化した値. なお, m ap40 とは, 観測結果として得られGroup Lasso の更新式に従い, w の更新を行う. 特に, 4 行

た画像上の対象天体を中心とした直径４秒角の円から得ら目の ∇f (w) は, 損失が小さくなる勾配を示している. また, 6

れる光量であり, 図 1 の赤いデータ点に対応している. こ行目の max(0, 1 − ηλ/||ŵGj ||)ŵGj で 0 が採択されない値まの値は, 直径 4 秒角の円内に存在する光量を全てを積算し

で wGj が大きくなった時, その wGj に対応したグループをモて値を算出しているため, 近隣の影響により, 対象天体が持デリングに重要な特徴量群とみなす. 上記の更新を収束まで繰

つ本来の光度とは離れた値が観測されてしまっている天体

り返し, 収束後の w についての決定係数 reg score を算出する.

データも存在するという問題点が存在する.

ここで得られた w 及び reg score は以下で説明する分類タスクに使用する特徴量の一部として扱う.

z apertureﬂux 40 mag 観測結果として得られた複数の画像を合成することで得られる Coadd 画像上の対象天体を中心とした直径４秒角の円から得られる光量. なお, Coadd
画像は, 対象天体が映っている全画像を画像的に足し合わ

1：{j, m + j} ∈ G は, それぞれ cos(2πti νj ), sin(2πti νj ) と対応している.

せることで得られる. 平均 m ap40 と比較し近傍天体の明diff_ap40の絶対値が⼤きい程, ノイズの影響を受けている
可能性が⾼いAlgorithm 2 VariableStarsExtraction
Input: i) データセット: Yii) 周波数候補数: miii) ナイキスト周波数 (表現可能な最大周波数): fn

23iv) ドメイン特徴量リスト (各天体が持つ固有の天体情報): Fd
Output: 変動天体リスト: L

平均m_ap40

22

1: // 変数設定
2: for i = 1 to n do

21

3:

for j = 1 to m do// ti = i,4:20

X[i][j] ← cos(2πti νj );

5:

19

6:

νj = j fmn
X[i][j + m] ← sin(2πti νj )end for

7: end for

18

8: for j = 1 to m do

17

10: end for

9:

-8

-6

-4

-2

0

diff_ap40
図 4: diﬀ ap40 と平均 m ap40 の関係

G[j] ← [j, j + m]

11: // 周期解析
12: for i = 1 to N do
13:

∗w, reg score ← Algorithm1(X, Y[i], G)

14:F [i] ← [∗w, reg score, Fd [i]]15: end forるさの影響を軽減する等, ある程度のノイズの除去が成さ

16: // 二値分類

れているため, 天体本来の明るさに近い値であることが期

17: F → Xtrain , Xtest , ytrain , ytest

待される.

18: lgb = LightGBM.train(Xtrain , ytrain )19: L = lgb.predict(Xtest )diﬀ ap40 平均 m ap40 と z apertureﬂux 40 mag の差. す

20: return Lなわち, 天体が影響を受けているノイズの大きさに相当する
特徴量. 具体的には, 平均 m ap40, z apertureﬂux 40 mag
が ど の 程 度 異 な る 値 で あ る か を 示 す 値 で あ る た め,z apertureﬂux 40 mag が天体本来の光度に近いと仮定
すると, diﬀ ap40 の絶対値が大きい程, ノイズとして考えられる近傍の天体の明るさの影響が大きいことを表す. 周
期解析では, 天体の光度（m ap40）の変動情報を用いているため, 平均 m ap40 に影響を与えているノイズの大きさ
は変動天体の識別に影響していると考えられる. そのため,天体が影響を受けているノイズの大きさを表す本特徴量を用いている.

data size 1 つの天体データに含まれる観測データ点の数. ス
パースモデリングを利用しているものの, 観測回数が少な
い一定数の天体データについては過学習を引き起こしてしまう. 従って, 周期解析結果が過学習を引き起こしている
蓋然性に相当する数値だと考え, 特徴量として利用することにした.

する. 2-11 行目では, 必要な変数を作成する. 具体的には, 2-8行目では説明変数 X を, 9-11 行目ではグループセット G を作成している. その後, 13-16 行目では, 周期解析を実行し, それらの結果及び Fd から二値分類の際に必要な特徴量リスト F を作成する. 最後に 18-20 行目では, 本章で説明した二値分類を実行している. 特に 18 行目で実行している訓練データ及びテ
ストデータへの分割については, テストデータに教師ラベルが振り分けられていないデータセットを用いることで, 変動天体の予測に活用することが可能である.

4 実

験

本章では, 実際の天体データを利用した実験を通じて, 既存手法との比較を行うことにより, 提案手法がより高精度な変動天体検出が可能であることを示す. また, Group Lasso を活用した周期解析により正しい周期情報の算出を行えていること, 正
しいドメイン知識を特徴量エンジニアリングできていることを

以下では, 各特徴量間の関係性について説明を行う. 図 4 で明確化する. 加えて, 提案手法がノイズが混入したデータに対は, diﬀ ap40 と平均 m ap40 の関係を表している. この図から,しても正しい周期解析が実行できることを示す.

diﬀ ap40 の絶対値が極端に大きい値を持つ天体 (図の左側) の多くが暗い天体であることがわかる. 故に, 暗い天体の方が外部からの影響をより受けやすいという背景が読み取れる. このような背景を変動天体検出に反映させることでより高度な分類を実現する.

4. 1 実 験 設 定
4. 1. 1 データセット
本研究では, 実天体データセットとして PDR2 2 (Public Data

Release 2) [13] を利用した. このデータセットのカラムは, 天

3. 2. 2 アルゴリズム
提案モデル全体を表現した擬似コードを Algorithm2 に記載

2：https://hsc-release.mtk.nao.ac.jp/doc/

体識別子 objectid, 観測時刻 mjd, 天体の光度 m ap40, 変動天

表 5: 各評価指標でのモデル評価

体識別子 objclass により構成される. また PDR2 に属する天

モデル

AUC

accuracy

precision

recallF1

体データ全てに対して objclass が割り振られているわけではな

B.Sesar 手法

0.750

0.743

0.507

0.507

0.507

線形回帰 + 無 4

0.633

0.729

0.432

0.109

0.168

線形回帰 + 有

0.846

0.820

0.702

0.538

0.605

Lasso + 無

0.854

0.825

0.769

0.478

0.586

4. 2. 2 節では, 周期解析のノイズ及びスパースに対するロバス

Lasso + 有

0.937

0.897

0.815

0.780

0.795

ト性を検証するために人工データを扱った実験を行った. 具体

Group Lasso + 無

0.860

0.828

0.793

0.476

0.590

Group Lasso + 有

0.939

0.899

0.820

0.785

0.801

く, 既に割り振られているデータの総数は計 1,557 件である. これらは, 既存手法の一つである B.Sesar らの手法 [2] で上位を占めるデータであり, 変動天体は約 26.1%を占めている. また,的なデータの含有量については, 最小で 0.01% ほどである. な3

提案手法

お, 本稿におけるデータの含有量とは, 時系列全体に対して欠損していない観測データ点の割合を示す.

4. 1. 2 精 度 評 価れると推察する.

4. 2. 2 周期解析のノイズ及びスパースに対するロバスト性の検証

ベースラインとして B.Sesar らの手法を利用する. こちらは現状国立天文台の方で変動天体検出のために使用されている手

初めに定性的な評価を行う. 図 6 は, 用意した人工データ (a),法であるため, ベースラインとして妥当であると考えた. 但し,(b) から約 0.01% 5 ほどのデータ点を抽出し作成したスパースな

この手法は教師あり学習による手法でないため, 提案モデルとデータに対して線形回帰及び Group Lasso を利用し, 周期解析

同等の方法である k 分割交差検証を用いた精度評価ができなを行った結果を描画したものである. (c), (d) が Group Lasso

い. そのため, データセット全体を占める変動天体の割合が約により得られた結果であり, (e), (f) が線形回帰により得られた

26.1%であることから, 二値分類タスクによって得られる予測結果である. 使用データ量が極端に少ないため, 線形回帰を活用

結果を占める変動天体の割合も同等であると仮定すると, ベーし算出した回帰曲線は, 技術的課題点に挙げた通り過学習が生

スライン手法によってスコアリングした結果の上位 26.1%を変じている. また, 特徴量選択が行われていないため, 変動天体の動天体候補として考えることによって, 比較的公平な条件で精重要周期を表現できていないことも読み取れる. 一方, Group

度評価が可能であると考えた. また, 本稿で扱う 1,557 件のデー

Lasso により得られたパワースペクトルからは重要特徴量が一タは光度が相対的に高いレコードの集合であるため, ノイズの点に集約されていることが汲み取れる. 加えて, 正解データと影響は無視しても良い程度である. 故に, 主要因である標準偏Group Lasso により得られた回帰曲線が類似していることも明差によりスコアリングを行う.

白である.
続いて定量的な評価を行う. 図 7 は使用データ点の含有率と,4. 2 実 験 結 果

線形回帰及び Group Lasso を使用することで得られる回帰曲

4. 2. 1 実天体データへの適用

線と正解データの整合率の関係を表した折れ線グラフである.

表 5 は, objclass 割り振られている全データ計 1,557 件に対評価指標は, 時系列間の周期性の整合度合いを計測可能なコサ

して, 4. 1. 2 節で仮定した条件のもと, 既存手法及び提案手法

イン類似度 [14] を用いた. 他の評価値として, DTW(Dynamic

を適用することで得られた評価値 AUC, accuracy, precision,

Time Warping) [13] が存在するが, こちらは時系列間の周期でrecall, F1 をまとめたものである. 提案手法は, ベースライン手

はなく距離に基準を設けているため, 本実験にはそぐわない評法と比較し, 大幅な性能向上を示す結果が見られる. また, ス

価値であると考えた. 実験結果に着目すると, 線形回帰を利用しパースモデリング及びドメイン知識の有効性を示すために, そた周期解析では, データの含有率が 1%を下回ったあたりから,れらに関する対照実験も行った. 結果は, 線形回帰を使用して整合度が急激に低下し, 誤差も大きくなっている. 但し, 提案手いるモデルや, ドメイン知識を活用していないモデルと比較し,法では, データの含有量が実際の天体データの平均的な欠損値

スパースモデリングの一種である Group Lasso を使用し, かつ,割合である 0.01%程になったとしても, 整合度はほぼ 1 を保っ

ドメイン知識を活用したケースが最も高精度に変動天体を検出ている. また, 誤差も小さく安定している. これらの結果から,できていることを示すことができた. ここから, 特徴量選択を正提案手法の周期解析にロバスト性があるということが言える.

しいグループ単位で行い, 得られた周期情報をうまく活用でき4. 2. 3 教師ラベルが未定な天体データへの応用ていることがわかる. また, 周期解析を行う際, Lasso を使用し

図 8 は教師ラベルが未定なレコードが含まれた天体データ計たモデルと Group Lasso を使用したモデルの性能差がほとん

5,478 件に対して提案モデルを適用し, 得られた変動天体候補

ど見られない. こちらは, 今回はベースライン手法で得られたリストの上位 1,000 件を自然科学研究機構国立天文台の方々に結果の上位を占めているデータセットを使用しているため, 特調査していただいた結果を, 横軸が上位 k 件, 縦軸がその k 件

徴量のグループ化を行わずともほとんど望ましい結果が得られ内の変動天体の割合という形で描画したものである. 但し, 学たためであると考える. ベースライン手法で取得しきれなかっ習に用いたデータをテストデータにも含めた形式で予測を行った天体データも含めて同様の実験を行うと, より性能の差が現
4：有, 無とはそれぞれドメイン知識を活用したか否かを示している.
3：天体データの平均的な欠損値数に基づいている

5：天体データの平均的な欠損値数に基づく

時間 [⽇]

×10!"

パワー

等級 等級

周波数 [Hz]時間 [⽇]×10!"

周波数 [Hz](a) 正解データ(b) 正解パワースペクトル

等級

パワー

時間 [⽇]×10!"

周波数 [Hz](c) Group Lasso にて算出した回帰曲線(d) Group Lasso にて算出したパワースペクトル
パワー

等級 等級

時間 [⽇]

×10!"

時間 [⽇]

周波数 [Hz]時間 [⽇](e) 線形回帰にて算出した回帰曲線(f) 線形回帰にて算出したパワースペクトル

図 6: 周期解析結果の比較
等級(a) の人工データを元に周期解析を行い (c)-(f) を算出した. (c), (e) に含まれている赤点が (a) の内, 周期解析に使用したデータ計 30 点である.
また, (a) のグラフに描かれている灰色部は元のデータに付与したノイズである.

評価値
0.9906 vs 0.1342
等級

1.0

パワー

時間 [⽇]

97.8 %

100×10!"

周波数 [Hz]500

0.8
時間 [⽇]

変
動
天
体
割
合

400

変
動
天
体

80
70

]

%

60

件

]

53.6 %

[

0.2

90

[

コ
サ
イ 0.6
ン
類 0.4
似
度

323 件

0.0

50
1.00

0.10

0

200

400

600

800

1000

変動天体候補 top-k [件]17

図 7: データ含有量と正解データとの整合度の関係

図 8: 候補リスト上位 k 件を占める変動天体の割合

た. このような実験設定にした主な理由は, 訓練データから得データが十分でない環境であったため, 適切な割合で提案モデられるスコアという基準と予測対象の差, すなわち, 変動天体でルを活用することでより高度な検出が期待できる.

ある確率に関連した情報を得ることができるためである. 図より, 上位 323 件内の変動天体の割合は約 97.8%である. この区分の大半は学習に使用した変動天体データによって構成されているので, ここに属する教師ラベルが未定なデータに関しては高確率で変動天体であることが言える. また, 上位 1,000 件内の変動天体の割合は約 53.6%である. 但し, 予測対象全件の教師ラベルが明確化されているわけではないため, データセットに属する全変動天体の内, どの程度検出したかどうかが不明である点に注意していただきたい. ところで, 4. 1. 1 節で説明した通り, ベースライン手法によって取得した上位 1,557 件内の変動天体の割合は約 26.1%である. それを踏まえると, 上位 1,000件を確認するだけで 200 件ほど追加で変動天体の検出できたことがわかる. また, 一般的には訓練データとテストデータの割
合が 7:3 程度であるのに対して, 本実験では約 3:7 という訓練

200
100
0

0.01

データの含有率 [%]300

5 関連研究5. 1 変動天体検出
変動天体の検出に関する代表的な手法として, B.Sesar らによる手法 [2] が挙げられる. これは, 提案モデルのような機械学
習によるアプローチとは異なり, intrinsic variability σ, χ2 , 等級の測定値の平均, 標準偏差, 歪度などといった統計的数値を活用したスコアリングを行うことで変動天体の検出を行っている.
変動天体かどうかの閾値は主に σ である. しかし, 天体の持つ明るさが暗ければ暗いほどノイズの影響を受けやすく, 例え σが基準を満たしていたとしても, それが変動天体固有の特徴によるものか, ノイズの影響によるものかが判断つかない. このようなケースに対応するため, χ2 を用いて解決している.

0こちらの手法では変動天体の持つと言われている周期情報を活用できていないため, 与えられたデータセットを満足に活用できていない.
[5]

5. 2 周 期 解 析
周期解析とは, 正弦波や余弦波など周期性を持った関数を基底とし, 入力データを近似・分析を行う手法のことである. また, 今[6]回必要な情報は入力データに潜む重要周期である. CLEAN [15]は雑多なパワースペクトル情報から重要な周期長の抽出を可能[7]にする. しかし, 入力データから重要周期を直接計算できるわけではない. また, lomb-scargle periodogram [5, 6] では, 最小[8]二乗フィッティングを活用することによって特定の周期長に対応したフーリエ振幅を計算するピリオドグラム法の一種である.
サンプリング間隔が一定でない時系列に対しても適用可能であ[9]る反面, ノイズの影響が大きい場合や欠損値がある程度含まれた時系列からは望ましい結果が得られない. また, 複数の周期性が混在したような時系列に対して適用した場合, 複数の周期[10]性の内一つしか取得できない場合や, 存在する周期性が互いに[11]

影響を及ぼし合った結果存在しない周期性が誤検出される場合がある, といった問題点が挙げられる.

6 おわりに[12][13]

本稿では, 周期性に着目した変動天体の検出手法を提案した.
提案モデルは, 変動天体が持つとされている周期性に着目し, 特定の周波数を指定することなく, 加えて, スパースモデリングの恩恵により入力データに潜むノイズの影響も軽減し, 変動天体

[14]を検出することが可能である. 実験では, 実際の天体データに対しても既存手法と比べてより高精度な結果が得られることを示した. また, ノイズを加えたスパースな人工データに対しても,[15]

正確な周期解析を行うことが可能なことを明確にした. 今後の課題としては, 並列処理を導入した高速化の実現や他のスパー
スモデリングを活用したモデルの性能改善などが考えられる.

7 謝

[16]

artowicz, Grégory Jevardat de Fombelle, Isabelle LecoeurTaı̈bi, Lea Karbevska, et al. Gaia data release 3: All-skyclassiﬁcation of 12.4 million variable sources into 25 classes.
arXiv preprint arXiv:2211.17238, 2022.
Nicholas R Lomb. Least-squares frequency analysis ofunequally spaced data. Astrophysics and space science,39(2):447–462, 1976.
Jeﬀrey D Scargle. Studies in astronomical time series analysis. ii-statistical aspects of spectral analysis of unevenly
spaced data. AJ, 263:835–853, 1982.
Robert Tibshirani. Regression shrinkage and selection via
the lasso. Journal of the Royal Statistical Society: Series
B (Methodological), 58(1):267–288, 1996.
Ming Yuan and Yi Lin. Model selection and estimation in regression with grouped variables. Journal of the Royal Statistical Society: Series B (Statistical Methodology), 68(1):49–
67, 2006.
Guolin Ke, Qi Meng, Thomas Finley, Taifeng Wang, Wei
Chen, Weidong Ma, Qiwei Ye, and Tie-Yan Liu. Lightgbm:
A highly eﬃcient gradient boosting decision tree. Advancesin neural information processing systems, 30, 2017.
Jerome H Friedman. Greedy function approximation: agradient boosting machine. AS, pages 1189–1232, 2001.
Tianqi Chen and Carlos Guestrin. Xgboost: A scalable tree
boosting system. In Proceedings of the 22nd acm sigkdd
international conference on KDD, pages 785–794, 2016.
Taichi Kato and Makoto Uemura. Period analysis using
the least absolute shrinkage and selection operator (lasso).
PASJ, 64(6), 2012.
Hiroaki Aihara, Yusra AlSayyad, Makoto Ando, Robert
Armstrong, James Bosch, Eiichi Egami, Hisanori Furusawa,
Junko Furusawa, Andy Goulding, Yuichi Harikane, et al.
Second data release of the hyper suprime-cam subaru strategic program. PASJ, 71(6):114, 2019.
Anna Huang et al. Similarity measures for text document clustering.
In Proceedings of the sixth newzealand computer science research student conference(NZCSRSC2008), Christchurch, New Zealand, volume 4,
pages 9–56, 2008.
David H Roberts, Joseph Lehár, and John W Dreher. Time
series analysis with clean-part one-derivation of a spectrum.
AJ, 93:968, 1987.
Ernest K Ryu and Wotao Yin. Proximal-proximal-gradientmethod. arXiv preprint arXiv:1708.06908, 2017.

付

辞

この成果は，国立研究開発法人新エネルギー・産業技術総合

A 補題 1 の証明

開発機構（ＮＥＤＯ）の委託業務（JPNP16007) の結果得られ証明.

たものです.

うになる.
文

録

ELasso を w について微分を行った結果は (A.1) 式のよ

献[1] AC Becker, JJ Bochanski, SL Hawley, Ž Ivezić, AF Kowalski, B Sesar, and AA West. Periodic variability of low-massstars in sloan digital sky survey stripe 82. AJ, 731(1):17,
2011.
[2] Branimir Sesar, Željko Ivezić, Robert H Lupton, Mario
Jurić, James E Gunn, Gillian R Knapp, Nathan De Lee,
J Allyn Smith, Gajus Miknaitis, Huan Lin, et al. Exploring the variable sky with the sloan digital sky survey. AJ,
134(6):2236, 2007.
[3] Eran O Ofek, Maayane Soumagnac, Guy Nir, Avishay GalYam, Peter Nugent, Frank Masci, and Shri R Kulkarni. A
catalogue of over 10 million variable source candidates in
ztf data release 1. MNRAS, 499(4):5782–5790, 2020.
[4] Lorenzo Rimoldini, Berry Holl, Panagiotis Gavras, Marc
Audard, Joris De Ridder, Nami Mowlavi, Krzysztof Nien-∂ELasso1= X ⊤ (Xw − y) + λsign(wk )
∂w
n

(A.1)

よって, ELasso を wk について微分を行った結果は (A.1) 式で得られたベクトルの k 番目を考えれば良いため,
nm∑∂ELasso1∑=xik (wj xij − yj ) + λsign(wk ) (A.2)
∂wk
n i=1
j=1

ELasso を最小化する wk を計算するためには, ELasso が凸関数であるため, (A.2) 式が 0 になるケースを考えれば良い. このようなケースに対して, wk を計算すると, (A.3) 式のように

なる.

∑n
wk =

i=1 (yi −∑mj =k| wj xij − w0 )xik − nλsign(wk )∑n2i=1 xik(A.3)
これらの計算により ELasso を最小化する w を計算できた.
但し, (A.3) 式は −nλ < wk < nλ について wk を定義できな
い. すなわち, この範囲においては微分不可能であることを意味している. ここで劣微分の考えを適用すると, この範囲において wk = 0 と考えても良いため, 証明が完了した.

B 補題 2 の証明
証明.

近接勾配法 [16] を活用する. これは, 微分不可能点が含まれている関数に対して勾配法を適用できるように改良した手法である. 近接勾配法の更新は以下の規則に従い行われる.

xi+1 = proxηg (xi − η∇f (xi ))1||x − y||2 }2proxg (y) = argmin{g(x) +
x

(B.1)

Group Lasso の損失関数である EGroupLasso を (B.1) に従い変形すると, 以下のように書き下せる.

wi+1 = argmin{h(w)}w

h(w) = ∇f (wi )⊤ (w − wi ) +f (w) = ||y − Xw||2 , g(w) =1||w − wi ||2 + λg(wi )2ηJ∑||wGj ||

j=1

すなわち, h(w) を最小化する w の内, w = 0 となる ŵGjの範囲が ||ŵGj || < ηλ であることを示すことができれば良い.
続いて, あるグループ Gj に対して, wGj に関する h(w) の微分を考えると, 以下の式が成立する.


w Gj

wGj = ŵi Gj − ηλ||wGj ||
||ŵi Gj || <= ηλ(wGj =| 0)(B.2)(wGj = 0)(B.2) 式を解くことで, h(w) を最小化する wGj を求めることができる. 変形することで以下の式が得られる.

wGj = max(0, 1 −ηλ)ŵi Gj||ŵi Gj ||(B.3)(B.3) より, 1 − ||ŵηλ< 0 の場合, wGj = 0 となり, その||iGj

範囲が ||ŵGj || < ηλ であるため, 証明が完了した.