WO2022079904A1 - パラメータ推定装置、パラメータ推定システム、パラメータ推定方法、及びプログラム - Google Patents

パラメータ推定装置、パラメータ推定システム、パラメータ推定方法、及びプログラム Download PDF

Info

Publication number
WO2022079904A1
WO2022079904A1 PCT/JP2020/039119 JP2020039119W WO2022079904A1 WO 2022079904 A1 WO2022079904 A1 WO 2022079904A1 JP 2020039119 W JP2020039119 W JP 2020039119W WO 2022079904 A1 WO2022079904 A1 WO 2022079904A1
Authority
WO
WIPO (PCT)
Prior art keywords
calculation
parameter estimation
time
time point
time points
Prior art date
Application number
PCT/JP2020/039119
Other languages
English (en)
French (fr)
Inventor
気吹 三品
浩気 濱田
大 五十嵐
亮 菊池
Original Assignee
日本電信電話株式会社
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by 日本電信電話株式会社 filed Critical 日本電信電話株式会社
Priority to EP20957733.7A priority Critical patent/EP4231272A4/en
Priority to JP2022556813A priority patent/JP7456514B2/ja
Priority to US18/245,195 priority patent/US20230367846A1/en
Priority to AU2020472128A priority patent/AU2020472128B2/en
Priority to CN202080106085.5A priority patent/CN116324935A/zh
Priority to PCT/JP2020/039119 priority patent/WO2022079904A1/ja
Publication of WO2022079904A1 publication Critical patent/WO2022079904A1/ja

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G06F17/16—Matrix or vector computation, e.g. matrix-matrix or matrix-vector multiplication, matrix factorization
    • H—ELECTRICITY
    • H04—ELECTRIC COMMUNICATION TECHNIQUE
    • H04L—TRANSMISSION OF DIGITAL INFORMATION, e.g. TELEGRAPHIC COMMUNICATION
    • H04L9/00—Cryptographic mechanisms or cryptographic arrangements for secret or secure communications; Network security protocols
    • H04L9/08—Key distribution or management, e.g. generation, sharing or updating, of cryptographic keys or passwords
    • H04L9/0816—Key establishment, i.e. cryptographic processes or cryptographic protocols whereby a shared secret becomes available to two or more parties, for subsequent use
    • H04L9/085—Secret sharing or secret splitting, e.g. threshold schemes
    • H—ELECTRICITY
    • H04—ELECTRIC COMMUNICATION TECHNIQUE
    • H04L—TRANSMISSION OF DIGITAL INFORMATION, e.g. TELEGRAPHIC COMMUNICATION
    • H04L2209/00—Additional information or applications relating to cryptographic mechanisms or cryptographic arrangements for secret or secure communication H04L9/00
    • H04L2209/46—Secure multiparty computation, e.g. millionaire problem

Definitions

  • the present invention relates to a technique for realizing parameter estimation of a cox proportional hazard model while maintaining confidentiality in a secret calculation technique.
  • Cox proportional hazards regression which is a regression analysis using the proportional hazard model of cox, is an analysis method often used in survival time analysis (Non-Patent Document 1). Cox proportional hazard regression can be performed in plain text in a package of commercially available statistical analysis software or software language.
  • a method called secret calculation is known as a method of obtaining a specific calculation result without restoring the encrypted numerical value.
  • encryption is performed by distributing pieces of numerical values to multiple secret computing devices, and multiple secret computing devices perform cooperative calculation to add / subtract, constant addition, multiplication, and constant multiplication without restoring the numerical values.
  • Logical operations negative, logical product, logical sum, exclusive logical sum
  • data format conversion integer and binary number
  • the calculation of parameter estimation of the cox proportional hazard model includes many processes such as exponent, division, and group-by sum, which have high processing costs in secret calculation, so that it is difficult to calculate efficiently by secret calculation. be.
  • the present invention has been made in view of the above points, and an object of the present invention is to provide a technique for efficiently estimating parameters of a cox proportional hazard model without decoding a value at a time point.
  • a parameter estimation device that performs parameter estimation of the cox proportional hazard model by secret calculation.
  • a data storage unit that stores a database that has a record including the time when an event is observed, the feature amount of the observation target at that time, and the state of the observation target at that time for each observation target.
  • a substitution table and a flag indicating the boundary between time points are generated, and by using the substitution table and the flags, the value at the time point is obtained.
  • a calculation unit that aggregates the feature quantities at each time point and estimates the parameters based on the aggregated results, while keeping the information secret.
  • a parameter estimation device including an output unit that outputs the parameters estimated by the calculation unit is provided.
  • a technique for efficiently estimating the parameters of the cox proportional hazard model without decoding the value at the time point is provided.
  • FIG. 1 shows a configuration diagram of the parameter estimation device 100 according to the present embodiment.
  • the parameter estimation device 100 in the present embodiment has an input unit 110, a calculation unit 120, an output unit 130, and a data storage unit 140.
  • the parameter estimation device 100 may be configured by one device (computer) or may be configured as a system composed of a plurality of computers. This system may be called a parameter estimation system.
  • the arithmetic unit 120 and the data storage unit 140 may be separate servers.
  • Data concealed from the data obtained by observation is input to the input unit 110 of the parameter estimation device 100.
  • the data handled by the parameter estimation device 100 is concealed data, and the calculation is performed by secret calculation.
  • the input data is stored in the data storage unit 140 as a database.
  • the arithmetic unit 120 estimates the parameters of the cox proportional hazard model by performing the processing described later on the data such as scalars, vectors, and matrices read from the database of the data storage unit 140.
  • the output unit 130 outputs the parameters estimated by the calculation unit 120.
  • the parameters calculated by the calculation unit 120 may be stored in the data storage unit 140 and output from the output unit 130 according to the access from the outside. The details of the processing contents in the calculation unit 120 will be described later.
  • the parameter estimation device 100 in the present embodiment can be realized by, for example, causing a computer to execute a program describing the processing contents described in the present embodiment.
  • the "computer” may be a physical machine or a virtual machine on the cloud.
  • the "hardware” described here is virtual hardware.
  • the above program can be recorded on a computer-readable recording medium (portable memory, etc.), saved, and distributed. It is also possible to provide the above program through a network such as the Internet or e-mail.
  • FIG. 2 is a diagram showing an example of the hardware configuration of the above computer.
  • the computer of FIG. 2 has a drive device 1000, an auxiliary storage device 1002, a memory device 1003, a CPU 1004, an interface device 1005, a display device 1006, an input device 1007, an output device 1008, and the like, which are connected to each other by a bus B, respectively.
  • the program that realizes the processing on the computer is provided by, for example, a recording medium 1001 such as a CD-ROM or a memory card.
  • a recording medium 1001 such as a CD-ROM or a memory card.
  • the program is installed in the auxiliary storage device 1002 from the recording medium 1001 via the drive device 1000.
  • the program does not necessarily have to be installed from the recording medium 1001, and may be downloaded from another computer via the network.
  • the auxiliary storage device 1002 stores the installed program and also stores necessary files, data, and the like.
  • the memory device 1003 reads and stores the program from the auxiliary storage device 1002 when the program is instructed to start.
  • the CPU 1004 realizes the function related to the parameter estimation device 100 according to the program stored in the memory device 1003.
  • the interface device 1005 is used as an interface for connecting to a network.
  • the display device 1006 displays a GUI (Graphical User Interface) or the like by a program.
  • the input device 1007 is composed of a keyboard, a mouse, buttons, a touch panel, and the like, and is used for inputting various operation instructions.
  • the output device 1008 outputs the calculation result.
  • a: b
  • the vector symbol ( ⁇ ) to be placed above the beginning of the character is described before the character, such as “ ⁇ a”.
  • the tensor on the third floor is described in italicized characters.
  • the character of the tensor on the third floor is represented by adding a tensor to the upper left of the character, for example, tensor ZZ'.
  • the cox proportional hazard model is a model represented by the equation (1) (Non-Patent Document 1).
  • t, ⁇ , and z represent time, weight, and feature, respectively, and ⁇ 0 (t) and exp ( ⁇ ⁇ T ⁇ z) are the baseline hazard function and relative risk function (hazard), respectively. Is called.
  • the weight parameter ⁇ is estimated. The weight is calculated by performing maximum likelihood estimation using the partial likelihood shown in Eq. (2).
  • D is the number of time points at which death was observed
  • ⁇ zi represents the characteristic amount of the patient who died at time point i .
  • death is used as a target event, but this is an example.
  • the event may be for the purpose of falling, inability to walk, onset of illness, hospitalization, or the like. Censoring may be interpreted as one of the events.
  • R i in equation (2) is a set of patients who have not been censored or died until just before time point i, and is called a risk set.
  • "censoring" means that observation becomes impossible and it becomes unknown whether or not death has occurred after that.
  • the partial likelihood function of Eq. (2) is calculated by calculating (hazard of deceased patient) / (sum of hazards of risk set) for each time point and multiplying by all time points. Since this partial likelihood is based on the assumption that multiple censorships and deaths do not occur at the same time (no tie data), the actual data that often has tie data is the Breslow shown in equation (3). The method is often used.
  • the Newton method or the like is generally used as a method for obtaining the maximum likelihood estimator of ⁇ .
  • the Newton method is used.
  • Newton's method after the equation (3) is transformed into a log-likelihood function, the calculation is performed using the first-order derivative (gradient) and the second-order derivative (hesian) of the log-likelihood function.
  • the log-likelihood function l ( ⁇ ⁇ ), its first derivative U ( ⁇ ⁇ ), and its second derivative I ( ⁇ ⁇ ) are shown in Eqs. (4), (5), and (6), respectively.
  • Equation (7) converges after about 5 iterations.
  • a value in which a certain value a is concealed by encryption, secret sharing, or the like is called a ciphertext or a concealed value of a, and is described as [a].
  • [a] refers to a set of secret sharing fragments possessed by each secret computing device.
  • the parentheses "[", "]” indicating that they are ciphertexts are slightly different from the format of the parentheses in the mathematical formulas inserted in the drawings and the specification, but in the specification text, for convenience of description, "[ You are using ",”] ".
  • the process of dividing the ciphertext [a] by the plaintext b is written in a notation such as [a] / b.
  • the notation is similarly [ ⁇ a] + [ ⁇ b], [A] + [B].
  • prefix sub Vector [ ⁇ a]: ([a 1 ], [a 2 ], ..., [an]) and scalar [b] to ([b], [b]-[a 1 ], [b] -([ A 1 ] + [a 2 ]), ..., [b]- ⁇ [ ⁇ a]) is described as prefixSub ([ ⁇ a], [b]). ..
  • Group-by common Group-by common is a process for generating intermediate data that can be commonly used in various Group-by operations such as Group-by sum and Group-by count.
  • the intermediate data includes a substitution table [ ⁇ ⁇ ] and a flag [ ⁇ e] indicating whether or not the value of the key is a boundary, and by using these, various Group-by operations using the same key can be performed efficiently. be able to.
  • substitution table ⁇ ⁇ is a vector indicating which element in the vector input to the Group-by common should be moved to sort the elements of the vector.
  • Flag ⁇ e is different by putting 0 in the position of the element if it is the same as the value under the element for each element of the sorted vector of the vector input to Group-by common. In the case, it is a vector containing 1.
  • FIG. 4 shows an example of inputs and flags. As shown in FIG. 4, the last value of the flag cannot be compared with the value below, so 1 is entered.
  • the processing is performed for each column.
  • the size of the output becomes smaller than the size of the input, but in the present embodiment, the sizes of the input and the output are the same, and the unnecessary portion is padded with 0 at the end. This makes it possible to conceal the number of point-in-time key attributes.
  • the result of Group-by sum is described as "a vector whose length is the number of time points" or "a matrix of the number of time points x n", but it is not actually necessary. Is padded with 0, and it deals with "a vector whose length is the number of records" and "a matrix of the number of records x n".
  • the calculation unit 120 of the parameter estimation device 100 reads the ciphertext data stored in the data storage unit 140 (database), and calculates the above-mentioned equations (5), (6), and (7) by secret calculation. Then, the parameters of the cox proportional hazard regression are estimated. In the following, first, the characteristic operation will be described.
  • groupByCommon and groupBySum aggregation is performed for each time point while keeping the time point value secret, and unnecessary parts are padded at 0, so information on the number of time points is not leaked.
  • the length is treated as a vector of time points and calculated collectively, and when the value at each time point is a vector with length n, it is treated as a matrix of time points ⁇ n. Handle and calculate all together. If the value at each time point is an n ⁇ n matrix, it is treated as a third-order tensor with the number of time points ⁇ n ⁇ n and calculated collectively.
  • the parameter estimation device 100 minimizes costly processing such as exp, division, and group-by sum in cox proportional hazard regression, and calculates efficiently.
  • exp is required for 7 times ⁇ number of time points and division is required for 3 times ⁇ number of time points per iteration of Newton's method. It was kept to a minimum as shown below.
  • the group-by common is utilized to efficiently calculate the secret calculation cox proportional hazard regression.
  • groupBySum only the aggregation using the flag [ ⁇ e] indicating the boundary is performed.
  • the parameter estimation device 100 ⁇ Detailed processing details> Next, the detailed processing contents executed by the parameter estimation device 100 will be described.
  • the concealed observation data is stored as a database in the data storage unit 140, and the arithmetic unit 120 processes the data by secret calculation to set the parameters of cox proportional hazard regression. presume.
  • the processing operation the above-mentioned characteristic operation is performed.
  • FIG. 5 shows an image of the database to be processed by the calculation unit 120.
  • FIG. 5 shows the data in plain text for convenience of explanation, and also shows the sorted state in ascending order at the time point.
  • the patient, the feature amount of the patient, the time point, and the state are recorded in the database.
  • the patient, the feature amount of the patient, the time point, and the state are recorded in the database.
  • the number of time points D the number of patients m.
  • the calculation unit 120 By reading the data from the database, the calculation unit 120 holds the feature amount for all patients as the matrix Z of m ⁇ n, the time point as the time vector ⁇ t, and the state for all the patients as the state vector ⁇ c.
  • the time vector ⁇ t is first sorted to create a substitution table ⁇ ⁇ , and by reusing it, the feature quantity Z and the state vector ⁇ c ⁇ t are used as keys. Sorting only needs to be sorted based on the substitution table ⁇ ⁇ . That is, the cost is lower than that of normal sorting.
  • the calculation unit 120 of the parameter estimation device 100 performs parameter estimation on the data in the above database according to the procedure of the algorithm shown in FIGS. 6 and 7.
  • Line numbers are added to FIGS. 6 and 7 for the sake of explanation. In the following, the line number of the processing part will be regarded as a step number for explanation.
  • step 3 of the algorithm 1 of FIG. 6 the arithmetic unit 120 initializes the nth-order weight vector [ ⁇ ⁇ ] with [0].
  • step 4 by Group-by common, ⁇ t is sorted to create a substitution table ⁇ ⁇ , and a flag ⁇ e is created.
  • steps 5 and 6 the arithmetic unit 120 sorts [Z] and [ ⁇ c], respectively, to create [Z'] and [ ⁇ c'].
  • step 8 the calculation unit 120 creates [ Z'dead ] from [Z'] in which the features other than the features of the dead cases are set to 0, and in step 9, from the sum of the features of the dead cases at each time point by groupBySum. [S] is created.
  • step 11 the number of deaths [ ⁇ d] at each time point is calculated.
  • step 13 [ tensor ZZ'], which is an m ⁇ n ⁇ n tensor, is created, and [ ⁇ ⁇ ] is updated in steps 15 to 17.
  • step 16 The processing of calcGH in step 16 will be described with reference to FIG. 7.
  • the arithmetic unit 120 corresponds to [W'] corresponding to ⁇ zexp ( ⁇ ⁇ T ⁇ z) and [ tensor X'] corresponding to ⁇ z ⁇ z T exp ( ⁇ ⁇ T ⁇ z). Etc. are calculated.
  • the [ t ⁇ v'] corresponding to the exp ( ⁇ ⁇ T ⁇ z j ) calculated in step 4 is reused in the subsequent calculations.
  • the arithmetic unit 120 calculates [ ⁇ v psub ] corresponding to ⁇ j ⁇ Ri exp ( ⁇ ⁇ T ⁇ z j ) at each time point.
  • Each element of the vector [ ⁇ v psub ] of the time point length is a scalar value ⁇ j ⁇ Ri exp ( ⁇ ⁇ T ⁇ z j ).
  • the calculation here is not an iterative process for each number of time points, but a calculation that processes all records at once. In the calculation of [W psub ] and [ tensor X psub ] described below, all the records are processed together in the same manner.
  • the arithmetic unit 120 calculates [W psub ] corresponding to ⁇ j ⁇ Ri ⁇ z j exp ( ⁇ ⁇ T ⁇ z j ) at each time point.
  • the arithmetic unit 120 calculates [ tensor X psub ] corresponding to ⁇ j ⁇ Ri ⁇ z j ⁇ z j T exp ( ⁇ ⁇ T ⁇ z j ) at each time point.
  • step 20 the arithmetic unit 120 calculates [ ⁇ y] corresponding to the reciprocal of ⁇ j ⁇ Ri exp ( ⁇ ⁇ T ⁇ z j ). This is the only part of the reciprocal calculation.
  • the calculation unit 120 calculates the equation (5) which is a gradient.
  • [ ⁇ y] and [ ⁇ d] are vectors whose length is the number of time points
  • [G] and [W] and [S] are matrices of the number of time points ⁇ the number of features
  • the calculation result [ t ⁇ g] ] Is a horizontal vector whose length is the number of features. The sum of all time points is calculated by sum in step 25.
  • the arithmetic unit 120 calculates the equation (6) which is a hesian.
  • the [G tmp ] [G tmp ] T in step 29 corresponds to the above-mentioned AAT .
  • the technique according to the present embodiment by processing all the data at once without performing the iterative processing as in the conventional plaintext processing, the number of time points can be calculated without decoding.
  • processes with high processing costs in secret calculations such as division, exponent, and group-by sum are parallelized, and cox proportional hazard regression is performed efficiently while concealing data and the number of time points. You will be able to estimate parameters.
  • a parameter estimation device that executes parameter estimation of the cox proportional hazard model by secret calculation.
  • a data storage unit that stores a database that has a record including the time when an event is observed, the feature amount of the observation target at that time, and the state of the observation target at that time for each observation target.
  • a calculation unit that aggregates the feature quantities at each time point and estimates the parameters based on the aggregated results, while keeping the information secret.
  • a parameter estimation device including an output unit that outputs parameters estimated by the calculation unit.
  • the first calculation unit executes the calculation of a plurality of exps in the calculation formula used in the iterative calculation for parameter estimation by the calculation of exp once per iteration and the calculation using the calculation result. Parameter estimator described in section.
  • the first calculation unit executes the calculation of a plurality of inverse numbers in the calculation formula used in the iterative calculation for parameter estimation by the calculation of the inverse number once per iteration and the calculation using the calculation result. Item or the parameter estimation device according to the second item.
  • the calculation unit performs the calculation for each time point in the calculation formula used in the iterative calculation for parameter estimation collectively for all time points using a vector, a matrix, or a tensor.
  • the parameter estimation device according to any one of the above items.
  • a parameter estimation system that executes parameter estimation of the cox proportional hazard model by secret calculation.
  • a data storage unit that stores a database that has a record including the time when an event is observed, the feature amount of the observation target at that time, and the state of the observation target at that time for each observation target.
  • a substitution table and a flag indicating the boundary between time points are generated, and by using the substitution table and the flags, the value at the time point is obtained.
  • a calculation unit that aggregates the feature quantities at each time point and estimates the parameters based on the aggregated results, while keeping the information secret.
  • a parameter estimation system including an output unit that outputs parameters estimated by the calculation unit. (Section 6) A parameter estimation method performed by a parameter estimation device that performs parameter estimation of the cox proportional hazard model by secret calculation.
  • a vector consisting of time points is read from a database that has records including the time when an event is observed, the feature amount of the observation target at that time point, and the state of the observation target at that time point for each observation target, and the vector is sorted. By doing so, a substitution table and a flag indicating the boundary between time points are generated, and by using the substitution table and the flag, the value of the time point is kept secret and the feature amount is aggregated for each time point. , The calculation step for estimating the parameters based on the aggregated result, A parameter estimation method including an output step for outputting the parameters estimated by the calculation step. (Section 7) A program for making a computer function as each part in the parameter estimation device according to any one of the items 1 to 4.
  • Parameter estimation device 110 Input unit 120 Calculation unit 130 Output unit 140 Storage unit 1000 Drive device 1001 Recording medium 1002 Auxiliary storage device 1003 Memory device 1004 CPU 1005 Interface device 1006 Display device 1007 Input device 1008 Output device

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Mathematical Physics (AREA)
  • Data Mining & Analysis (AREA)
  • Computational Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Mathematical Analysis (AREA)
  • Mathematical Optimization (AREA)
  • Pure & Applied Mathematics (AREA)
  • Software Systems (AREA)
  • Databases & Information Systems (AREA)
  • Algebra (AREA)
  • General Engineering & Computer Science (AREA)
  • Computing Systems (AREA)
  • Signal Processing (AREA)
  • Computer Networks & Wireless Communication (AREA)
  • Computer Security & Cryptography (AREA)
  • Complex Calculations (AREA)
  • Management, Administration, Business Operations System, And Electronic Commerce (AREA)
  • Electrophonic Musical Instruments (AREA)
  • Stored Programmes (AREA)
  • Communication Control (AREA)
  • Information Retrieval, Db Structures And Fs Structures Therefor (AREA)

Abstract

ｃｏｘ比例ハザードモデルのパラメータ推定を、秘密計算により実行するパラメータ推定装置において、イベントが観測された時点と、当該時点の観測対象の特徴量と、当該時点の観測対象の状態とを含むレコードを、観測対象毎に有するデータベースを格納するデータ格納部と、前記データベースから、時点からなるベクトルを読み出し、当該ベクトルをソートすることにより、置換表と、時点の境目を示すフラグとを生成し、前記置換表と、前記フラグとを用いることにより、時点の値を秘匿したまま、前記特徴量の時点毎の集計を行い、集計結果に基づいて前記パラメータ推定を行う演算部と、前記演算部により推定されたパラメータを出力する出力部とを備える。

Description

パラメータ推定装置、パラメータ推定システム、パラメータ推定方法、及びプログラム
　本発明は、秘密計算技術において、秘匿性を保ったままｃｏｘ比例ハザードモデルのパラメータ推定を実現する技術に関連するものである。
　ｃｏｘの比例ハザードモデルを用いた回帰分析であるｃｏｘ比例ハザード回帰は、生存時間解析でよく用いられる分析手法である（非特許文献１）。市販の統計解析ソフトウェアやソフトウェア言語のパッケージにおいて、平文でｃｏｘ比例ハザード回帰を行うことができる。
　また、暗号化された数値を復元すること無く特定の演算結果を得る方法として、秘密計算と呼ばれる方法が知られている。一例として、複数の秘密計算装置に数値の断片を分散させるという暗号化を行い、複数の秘密計算装置が協調計算を行うことにより、数値を復元すること無く、加減算、定数加算、乗算、定数倍、論理演算（否定、論理積、論理和、排他的論理和）、データ形式変換（整数と二進数）等の結果を複数の秘密計算装置に分散された状態として得ることができる。
D. R. Cox. Regression Models and Life-Tables. Journal of the Royal Statistical Society. Series B (Methodological), Vol. 34, No. 2. (1972), pp. 187-220.
　ｃｏｘ比例ハザードモデルのパラメータ推定を行う際、平文では死亡や打ち切りが発生した時点毎に計算を繰り返す。しかし、秘密計算でも同様の方法を用いる場合、秘匿しておくべき時点の値を復号する必要があるという課題がある。また、ｃｏｘ比例ハザードモデルのパラメータ推定の計算では、指数、除算、ｇｒｏｕｐ－ｂｙ　ｓｕｍといった秘密計算での処理コストが大きい処理が多く含まれるため、秘密計算で効率良く計算するのが難しいという課題もある。
　本発明は上記の点に鑑みてなされたものであり、時点の値を復号することなく、効率良くｃｏｘ比例ハザードモデルのパラメータ推定を行うための技術を提供することを目的とする。
　開示の技術によれば、ｃｏｘ比例ハザードモデルのパラメータ推定を、秘密計算により実行するパラメータ推定装置であって、
　イベントが観測された時点と、当該時点の観測対象の特徴量と、当該時点の観測対象の状態とを含むレコードを、観測対象毎に有するデータベースを格納するデータ格納部と、
　前記データベースから、時点からなるベクトルを読み出し、当該ベクトルをソートすることにより、置換表と、時点の境目を示すフラグとを生成し、前記置換表と、前記フラグとを用いることにより、時点の値を秘匿したまま、前記特徴量の時点毎の集計を行い、集計結果に基づいて前記パラメータ推定を行う演算部と、
　前記演算部により推定されたパラメータを出力する出力部と
　を備えるパラメータ推定装置が提供される。
　開示の技術によれば、時点の値を復号することなく、効率良くｃｏｘ比例ハザードモデルのパラメータ推定を行うための技術が提供される。
本発明の実施の形態におけるパラメータ推定装置の構成図である。 装置のハードウェア構成例を示す図である。 置換表を説明するための図である。 フラグを説明するための図である。 データベースの例を示す図である。 パラメータ推定装置の処理手順の例を説明するための図である。 パラメータ推定装置の処理手順の例を説明するための図である。
　以下、図面を参照して本発明の実施の形態（本実施の形態）を説明する。以下で説明する実施の形態は一例に過ぎず、本発明が適用される実施の形態は、以下の実施の形態に限られるわけではない。
　（装置構成例）
　図１に、本実施の形態におけるパラメータ推定装置１００の構成図を示す。図１に示すように、本実施の形態におけるパラメータ推定装置１００は、入力部１１０、演算部１２０、出力部１３０、及びデータ格納部１４０を有する。
　パラメータ推定装置１００は、１つの装置（コンピュータ）で構成されてもよいし、複数のコンピュータからなるシステムとして構成されてもよい。このシステムをパラメータ推定システムと呼んでもよい。パラメータ推定システムにおいて、例えば、演算部１２０とデータ格納部１４０が別々のサーバであってもよい。
　パラメータ推定装置１００の入力部１１０に、観測により得られたデータから秘匿化されたデータが入力される。特に断らなければ、パラメータ推定装置１００により扱われるデータは秘匿化されたデータであり、計算は秘密計算でなされるものとする。
　入力されたデータはデータベースとしてデータ格納部１４０に格納される。演算部１２０は、データ格納部１４０のデータベースから読み出したスカラー、ベクトル、行列等のデータに対して後述する処理を行うことで、ｃｏｘ比例ハザードモデルのパラメータ推定を行う。出力部１３０は、演算部１２０で推定されたパラメータを出力する。
　なお、演算部１２０で計算されたパラメータがデータ格納部１４０に格納され、外部からのアクセスに応じて、出力部１３０から出力されてもよい。演算部１２０での処理内容の詳細については後述する。
　（ハードウェア構成例）
　本実施の形態におけるパラメータ推定装置１００は、例えば、コンピュータに、本実施の形態で説明する処理内容を記述したプログラムを実行させることにより実現可能である。なお、この「コンピュータ」は、物理マシンであってもよいし、クラウド上の仮想マシンであってもよい。仮想マシンを使用する場合、ここで説明する「ハードウェア」は仮想的なハードウェアである。
　上記プログラムは、コンピュータが読み取り可能な記録媒体（可搬メモリ等）に記録して、保存したり、配布したりすることが可能である。また、上記プログラムをインターネットや電子メール等、ネットワークを通して提供することも可能である。
　図２は、上記コンピュータのハードウェア構成例を示す図である。図２のコンピュータは、それぞれバスＢで相互に接続されているドライブ装置１０００、補助記憶装置１００２、メモリ装置１００３、ＣＰＵ１００４、インタフェース装置１００５、表示装置１００６、入力装置１００７、出力装置１００８等を有する。
　当該コンピュータでの処理を実現するプログラムは、例えば、ＣＤ－ＲＯＭ又はメモリカード等の記録媒体１００１によって提供される。プログラムを記憶した記録媒体１００１がドライブ装置１０００にセットされると、プログラムが記録媒体１００１からドライブ装置１０００を介して補助記憶装置１００２にインストールされる。但し、プログラムのインストールは必ずしも記録媒体１００１より行う必要はなく、ネットワークを介して他のコンピュータよりダウンロードするようにしてもよい。補助記憶装置１００２は、インストールされたプログラムを格納すると共に、必要なファイルやデータ等を格納する。
　メモリ装置１００３は、プログラムの起動指示があった場合に、補助記憶装置１００２からプログラムを読み出して格納する。ＣＰＵ１００４は、メモリ装置１００３に格納されたプログラムに従って、パラメータ推定装置１００に係る機能を実現する。インタフェース装置１００５は、ネットワークに接続するためのインタフェースとして用いられる。表示装置１００６はプログラムによるＧＵＩ（Ｇｒａｐｈｉｃａｌ　Ｕｓｅｒ　Ｉｎｔｅｒｆａｃｅ）等を表示する。入力装置１００７はキーボード及びマウス、ボタン、又はタッチパネル等で構成され、様々な操作指示を入力させるために用いられる。出力装置１００８は演算結果を出力する。
　（準備）
　パラメータ推定装置１００の動作例を説明する前に、その準備として、記法、ｃｏｘ比例ハザード回帰、及び秘密計算等について説明する。
　　＜記法＞
　ａをｂで定義することをａ：＝ｂと記載し、ベクトルを→ａ：＝（ａ０，...，ａｎ－１）と記載し、特筆しない限りＡのような大文字は行列を表し、その転置行列はＡＴと記載する。なお、本明細書のテキストにおいて、記載の便宜上、文字の頭の上に置かれるべきベクトルの記号（→）を、「→ａ」のように、文字の前に記載している。また、図面に記載したアルゴリズムにおいて、３階のテンソルは斜体の文字で記載している。本明細書のテキストにおいては、記載の便宜上、３階のテンソルの文字を、例えばｔｅｎｓｏｒＺＺ´のように、その文字の左上にｔｅｎｓｏｒを付けて表す。
　加減乗算において入力がベクトル→ａもしくは行列Ａと、スカラーｂの場合、→ａ、Ａの全ての要素に対してｂとの演算を行うものとする。また、特に記載が無いベクトルは列ベクトルである。横ベクトルの場合はｔ→ａのように左上にｔを付けることで区別する。
　　＜ｃｏｘ比例ハザード回帰＞
　ｃｏｘ比例ハザードモデルは式（１）で表されるモデルである（非特許文献１）。
Figure JPOXMLDOC01-appb-M000001
 　式（１）において、ｔ、β、ｚはそれぞれ時間、重み、特徴量を表し、λ０（ｔ）、ｅｘｐ（→βＴ→ｚ）はそれぞれベースラインハザード関数、相対危険度関数（ハザード）と呼ばれる。ｃｏｘ比例ハザード回帰では重みのパラメータβを推定する。重みは、式（２）に示す部分尤度を用いて最尤推定を行うことで計算する。
Figure JPOXMLDOC01-appb-M000002
 　式（２）において、Ｄは死亡が観測された時点の数であり、→ｚｉは時点ｉに死亡した患者の特徴量を表す。なお、本実施の形態では、目的となるイベントとして、死亡を用いているが、これは例である。例えば、転倒、歩行不可能、疾患発症、入院等を目的のイベントとしてもよい。打ち切りをイベントの１つと解釈してもよい。
　式（２）におけるＲｉは、時点ｉの直前まで打ち切りも死亡も発生していない患者の集合であり、リスクセットと呼ばれる。なお、「打ち切り」とは、観測ができなくなり、それ以降、死亡の発生の有無が不明になることである。
　従って、式（２）の部分尤度関数は、時点毎に（死亡した患者のハザード）／（リスクセットのハザードの総和）を計算し、全時点分掛け合わせたものである。この部分尤度は同じ時刻に複数の打ち切りや死亡が発生していない（タイデータが無い）という仮定を置いているため、タイデータがあることの多い実データでは、式（３）に示すＢｒｅｓｌｏｗ法が良く用いられる。
Figure JPOXMLDOC01-appb-M000003
 　式（３）は、基本的には式（２）のｃｏｘの部分尤度と同じであるが、分母がｄｉ乗（ｄｉは時点ｉの死亡患者数）されている点と、分子に→ｚｉの代わりに、時点ｉの死亡患者の特徴量の総和ｓｉを用いる点が異なる。ｄｉ＝１の場合は式（２）と一致するため、式（３）は、式（２）を一般化した形と考えることができる。以下の説明では、式（３）を前提としている。
　→βの最尤推定量を求める方法としては、Ｎｅｗｔｏｎ法等が一般的である。本実施の形態ではＮｅｗｔｏｎ法を用いている。Ｎｅｗｔｏｎ法では、式（３）を対数尤度関数へと変形した後、その対数尤度関数の１階微分（勾配）と２階微分（ヘシアン）を用いて計算する。対数尤度関数ｌ（→β）、その１階微分Ｕ（→β）、及び２階微分Ｉ（→β）をそれぞれ式（４）、式（５）、式（６）に示す。
Figure JPOXMLDOC01-appb-M000004
Figure JPOXMLDOC01-appb-M000005
Figure JPOXMLDOC01-appb-M000006
 　Ｎｅｗｔｏｎ法では、式（５）と式（６）を用いて、以下の式（７）を反復して→βの最尤推定値を求める。式（７）は、およそ５回ほどの反復で収束する。
Figure JPOXMLDOC01-appb-M000007
 　＜秘密計算＞
　ある値ａを暗号化や秘密分散等により秘匿化した値をａの暗号文あるいは秘匿値と呼び、［ａ］と記述する。ａが秘密分散により秘匿化された場合は、［ａ］により各秘密計算装置が持つ秘密分散の断片の集合を参照する。なお、暗号文であることを示す括弧"［"、"］"について、図面及び明細書に挿入された数式での括弧の書式と少し異なるが、明細書テキストにおいては、記載の便宜上、"［"、"］"を使用している。
　次に、秘密計算における各種の演算を説明する。
　■四則演算
　２つの暗号文［ａ］、［ｂ］の加算、減算、乗算は、それぞれ暗号文［ａ＋ｂ］、［ａ‐ｂ］、［ａ×ｂ］を計算する処理である。これらの演算をそれぞれ、［ａ］＋［ｂ］、［ａ］‐［ｂ］、［ａ］×［ｂ］と記載する。
　また、暗号文［ａ］を平文ｂで割る処理は、［ａ］／ｂのような記法とする。入力がベクトルや行列で、要素毎にこれらの処理を行う場合も同様に［→ａ］＋［→ｂ］、［Ａ］＋［Ｂ］のような記法とする。
　加減乗算において入力が行列Ａと列ベクトル→ｂの場合は、行列の各列ベクトルに対して→ｂとの要素毎の演算を実施し、入力が行列Ａと行ベクトルｔ→ｂの場合は、行列の各行ベクトルに対してｔ→ｂとの要素毎の演算を実施するものとする。
　■総和
　ベクトル［→ａ］の要素の総和を求める処理をｓｕｍ（［→ａ］）と記述する。またｓｕｍ（［Ａ］）のようにｍ×ｎの行列がｓｕｍの入力の場合は、列方向の総和を計算し、長さｎの行ベクトル［ｔ→ｃ］を出力するものとする。
　■ｐｒｅｆｉｘ　ｓｕｂ
　ベクトル［→ａ］：＝（［ａ１］，［ａ２］，・・・，［ａｎ］）とスカラー［ｂ］から（［ｂ］，［ｂ］－［ａ１］，［ｂ］－（［ａ１］＋［ａ２］），・・・，［ｂ］－Σ［→ａ］）となるベクトルを計算することを、ｐｒｅｆｉｘＳｕｂ（［→ａ］，［ｂ］）と記載する。
　■逆数
　暗号文［ａ］の逆数１／［ａ］を計算することを、［ｃ］←ｒｅｃｉｐｒｏｃａｌ（［ａ］）のように記載する。入力がベクトルの場合も同じ記法とする。
　■指数
　暗号文［ａ］を入力とし、ネイピア数ｅの［ａ］乗を計算することを［ｃ］←ｅｘｐ（［ａ］）のように記載する。入力がベクトルの場合も同じ記法とする。
　■Ｇｒｏｕｐ－ｂｙ　ｃｏｍｍｏｎ
　Ｇｒｏｕｐ－ｂｙ　ｃｏｍｍｏｎは、Ｇｒｏｕｐ－ｂｙ　ｓｕｍやＧｒｏｕｐ－ｂｙ　ｃｏｕｎｔといった様々なＧｒｏｕｐ－ｂｙ演算で共通的に用いることができる中間データを生成する処理である。中間データは置換表［→π］と、キーの値の境目かどうかを表すフラグ［→ｅ］を含み、これらを使いまわすことで、同じキーを用いた様々なＧｒｏｕｐ－ｂｙ演算を効率良く行うことができる。
　図３、図４を参照して、Ｇｒｏｕｐ－ｂｙ　ｃｏｍｍｏｎの演算で得られる置換表→πと、フラグ→ｅの例を説明する。ここでは、説明の便宜上、平文の処理として説明する。
　置換表→πは、Ｇｒｏｕｐ－ｂｙ　ｃｏｍｍｏｎに入力されたベクトルにおける各要素が、何番目の要素に移動すれば、当該ベクトルの要素がソートされるかを表すベクトルである。
　図３において、ベクトル→ｋをＧｒｏｕｐ－ｂｙ　ｃｏｍｍｏｎに入力すると、当該ベクトルがソートされ、図３に示すとおりの置換表→πが得られる。例えば、置換表→πの最初の要素である３は、入力ベクトル→ｋの最初の要素が、３番目の要素になればソートできることを示している。
　フラグ→ｅは、Ｇｒｏｕｐ－ｂｙ　ｃｏｍｍｏｎに入力されたベクトルのソート済みのベクトルの各要素に対して、その要素の下の値と比較して同じならば、その要素の位置に０を入れ、異なる場合は１を入れたベクトルである。図４に、入力とフラグの例を示す。図４に示すように、フラグの一番最後については、下の値と比較できないので１を入れる。
　キーのベクトル［→ｋ］を入力してＧｒｏｕｐ－ｂｙ　ｃｏｍｍｏｎを行うことを式（８）のように記述する。置換表［→π］を用いてベクトル［→ａ］や行列［Ａ］（行数が［→π］の長さと等しい）をソートする処理を式（９）、式（１０）のように記述する。
　ソート済みのベクトル［→ａ´］や［Ａ´］と［→ｅ］を用いてＧｒｏｕｐ－ｂｙ　ｓｕｍを行うことを式（１１）（１２）のように記述する。［→ａ´］、［Ａ´］などのプライムは、ソート済みであることを表し、以降も同様の記法とする。
Figure JPOXMLDOC01-appb-M000008
Figure JPOXMLDOC01-appb-M000009
Figure JPOXMLDOC01-appb-M000010
Figure JPOXMLDOC01-appb-M000011
Figure JPOXMLDOC01-appb-M000012
 　ｓｏｒｔ、ｇｒｏｕｐＢｙＳｕｍの入力が行列の場合、処理は列ごとに行われる。また一般的にＧｒｏｕｐ－ｂｙ　ｓｕｍを行うと出力のサイズは入力のサイズ以下になるが、本実施の形態では入力と出力のサイズは同じであり、不要な分は末尾を０でパディングする。これによって時点キー属性の数を秘匿することができる。なお、後述する処理の説明において、便宜的にＧｒｏｕｐ－ｂｙ　ｓｕｍの結果を「長さが時点数のベクトル」や「時点数×ｎの行列」と記載しているが、実際には不要な部分を０でパディングした「長さがレコード数のベクトル」や「レコード数×ｎの行列」を扱っている。
　例えば、図４に示すフラグ→ｅを（０，１，０，０，１，１）と記載するとして、仮に→ａ´が（２，１，３，５，１，２）であるとすると、→ａ´と→ｅを入力とするＧｒｏｕｐ－ｂｙ　ｓｕｍは、（２＋１，３＋５＋１，２，０，０，０）＝（３，９，２，０，０，０）になる。
　（パラメータ推定装置１００の動作について）
　以下、パラメータ推定装置１００の動作例を説明する。パラメータ推定装置１００の演算部１２０は、データ格納部１４０（データベース）に格納された暗号文のデータを読み出し、秘密計算により、前述した式（５）、式（６）、式（７）を計算して、ｃｏｘ比例ハザード回帰のパラメータ推定を行う。以下では、まず、特徴的な動作について説明する。
　＜全レコードをまとめて処理する＞
　式（５）、式（６）をそのまま実装した場合、時点毎に反復計算していき、計算結果を順番に足していくような処理になる。秘密計算ｃｏｘ比例ハザード回帰を計算する場合には時点数を秘匿するため、前述のｇｒｏｕｐＢｙＣｏｍｍｏｎ、ｇｒｏｕｐＢｙＳｕｍを用いる。
　ｇｒｏｕｐＢｙＣｏｍｍｏｎ、ｇｒｏｕｐＢｙＳｕｍでは時点の値を秘匿したまま時点毎の集計を行い、不要な部分は０でパディングされるため、時点数に関する情報が漏れない。
　また、時点数毎の反復処理ではなく全レコードをまとめて処理することにより、秘密計算での処理コストが大きい演算の回数を減らすことができて処理効率も良い。つまり、例えば、時点１でΣｊ∈Ｒｉｅｘｐ（→βＴ→ｚｊ）を計算し、時点２でΣｊ∈Ｒｉｅｘｐ（→βＴ→ｚｊ）を計算し、...といった計算ではなく、全時点分のΣｊ∈Ｒｉｅｘｐ（→βＴ→ｚｊ）を一度に計算する。より具体的には、各時点の値がスカラーの場合は長さが時点数のベクトルとして扱ってまとめて計算し、各時点の値が長さｎのベクトルの場合は時点数×ｎの行列として扱って、まとめて計算する。また、各時点の値がｎ×ｎ行列の場合は時点数×ｎ×ｎの３階のテンソルとして扱ってまとめて計算する。
　＜処理コストの大きい演算の削減＞
　式（５）、式（６）にはｅｘｐや除算が多く含まれ、また、Σｊ∈Ｒｉの処理はＧｒｏｕｐ－ｂｙ　ｓｕｍであるため、秘密計算での計算コストが非常に大きい。
　本実施の形態に係るパラメータ推定装置１００は、ｃｏｘ比例ハザード回帰におけるｅｘｐ、除算、Ｇｒｏｕｐ－ｂｙ　ｓｕｍといったコストの大きい処理を最小限に抑え、効率良く計算する。単純に式（５）、式（６）の通りに計算した場合、Ｎｅｗｔｏｎ法の１反復あたりｅｘｐは７回×時点数、除算は３回×時点数分必要になるが、本実施の形態では下記の通り最小限に抑えた。
　・ｅｘｐの計算が１反復あたり１回
　・逆数の計算が１反復あたり１回
　また、Ｇｒｏｕｐ－ｂｙ　ｓｕｍについても、処理をｇｒｏｕｐＢｙＣｏｍｍｏｎ、ｇｒｏｕｐＢｙＳｕｍに分けることで、より効率的に処理を行うこととしている。これらの演算の削減についてより詳細に説明する。
　■ｅｘｐの削減
　式（５）、式（６）に示すとおり、ｅｘｐの引数が全て→βＴ→ｚｊであるため、一度計算したら後は使い回せば良い。これに加えて、前述の全レコードをまとめて処理することで処理が並列化され、１反復あたり１度で済むようになる。
　■除算の削減
　除算を逆数の計算＋乗算で行う場合、Σｊ∈Ｒｉｅｘｐ（→βＴ→ｚｊ）の逆数を、式（５）の第２項と式（６）の第１項で使い回せるため、除算２回ではなく逆数の計算１回＋乗算２回で済む。式（６）の第２項は前述の２つとは除数が異なるが、この項は除算をしなくても求めることができる。式（５）の第２項と式（６）の第２項を見比べてみると、式（５）の第２項のｄｉを除いた部分をＡとしたとき、式（６）の第２項はＡＡＴで表せるため、積のみで式（６）の第２項は計算できる。これに加えて、前述の全レコードをまとめて処理することで処理が並列化され、１反復あたり１度で済むようになる。
　■Ｇｒｏｕｐ－ｂｙ　ｓｕｍの削減
　式（５）、式（６）の通りに実装するとＧｒｏｕｐ－ｂｙｓｕｍを処理の中で何度も行うことになり、処理効率が低下する。そこで、本実施の形態では、全てキーが同一であることに着目し、最初にキーに対して一度だけ行ったＧｒｏｕｐ－ｂｙｃｏｍｍｏｎで得た［→ｅ］を使い回すこととしている。
　本実施の形態のパラメータ推定装置１００では、Ｇｒｏｕｐ－ｂｙ　ｃｏｍｍｏｎを活用し、秘密計算ｃｏｘ比例ハザード回帰を効率良く計算する。ｇｒｏｕｐＢｙＳｕｍの計算では、境目を表すフラグ［→ｅ］を用いた集計を行うのみである。
　＜詳細な処理内容について＞
　次に、パラメータ推定装置１００が実行する詳細な処理内容を説明する。ここでは、秘匿化された観測データがデータ格納部１４０においてデータベースとして格納されているものとし、演算部１２０は、秘密計算によりそのデータに対して処理を行うことで、ｃｏｘ比例ハザード回帰のパラメータを推定する。処理動作において、上述した特徴的な動作が行われる。
　演算部１２０による処理の対象となるデータベースのイメージを図５に示す。図５は、説明の便宜上、データが平文で示されており、また、時点の昇順にソート済の状態を示している。
　図５に示すように、データベースには、患者数ｍの患者（観測対象）毎に、ｎ個の特徴量と、観測の時点と、その時点での状態（死亡＝１、打ち切り＝０）が格納されている。ｍはレコード数でもあり、Ｄは時点数である。例えば、時刻ベクトル→ｔ＝（１，１，１，２）であれば、時点数は２である。
　図５の例では、１患者に対し、死亡があった場合に、その患者、その患者の特徴量、その時点、その状態がデータベースに記録される。また、その時点で打ち切りされている患者についても、その患者、その患者の特徴量、その時点、その状態がデータベースに記録される。
　同時点で、複数患者に対して死亡又は打ち切りが観測される場合がある。そのため、時点数Ｄ≦患者数ｍである。
　演算部１２０は、データベースからデータを読み出すことで、患者全員分の特徴量をｍ×ｎの行列Ｚとし、時点を時刻ベクトル→ｔとし、患者全員分の状態を状態ベクトル→ｃとして保持する。
　Ｚ、→ｔ、→ｃはいずれも初期状態においてはソートされていない。前述したように、Ｇｒｏｕｐ－ｂｙ　ｃｏｍｍｏｎにより、最初に時刻ベクトル→ｔをソートして置換表→πを作り、それを使い回すことで、特徴量Ｚや状態ベクトル→ｃの→ｔをキーとするソートは、置換表→πに基づいた並べ替えのみで済む。すなわち、通常のソートを行うより低コストで済む。
　＜処理手順＞
　パラメータ推定装置１００の演算部１２０は、上記のデータベースのデータに対して、図６、図７に示すアルゴリズムの手順に従って、パラメータ推定を行う。図６、図７には、説明のために行番号を付している。以下では、処理部分の行番号をステップ番号とみなして説明を行う。
　図６のアルゴリズム１のステップ３において、演算部１２０は、ｎ次の重みベクトルである［→β］を［０］で初期化する。ステップ４において、Ｇｒｏｕｐ－ｂｙ　ｃｏｍｍｏｎにより、→ｔをソートして置換表→πを作成するとともに、フラグ→ｅを作成する。
　ステップ５、６において、演算部１２０は、［Ｚ］、［→ｃ］をそれぞれソートして、［Ｚ´］、［→ｃ´］を作成する。
　演算部１２０は、ステップ８において、［Ｚ´］から死亡例の特徴量以外を０とした［Ｚ´ｄｅａｄ］を作成し、ステップ９において、ｇｒｏｕｐＢｙＳｕｍにより時点毎の死亡例の特徴量の総和からなる［Ｓ］を作成する。ステップ１１において、時点毎の死亡数［→ｄ］を算出する。ステップ１３において、ｍ×ｎ×ｎテンソルである［ｔｅｎｓｏｒＺＺ´］を作成し、ステップ１５～１７において［→β］の更新を行う。
　ステップ１６のｃａｌｃＧＨの処理については、図７を参照して説明する。ステップ３～６において、演算部１２０は、→ｚｅｘｐ（→βＴ→ｚ）に対応する［Ｗ´］と、→ｚ→ｚＴｅｘｐ（→βＴ→ｚ）に対応する［ｔｅｎｓｏｒＸ´］等を算出する。ステップ４で計算したｅｘｐ（→βＴ→ｚｊ）に対応する［ｔ→ｖ´］が以降の計算で使い回されている。
　ステップ８～１０において、演算部１２０は、各時点のΣｊ∈Ｒｉｅｘｐ（→βＴ→ｚｊ）に対応する［→ｖｐｓｕｂ］を算出する。時点数長のベクトルである［→ｖｐｓｕｂ］の各要素が、スカラー値であるΣｊ∈Ｒｉｅｘｐ（→βＴ→ｚｊ）になっている。つまり、ここでの計算は、時点数毎の反復処理でなく、全レコードをまとめて処理する計算である。下記で説明する［Ｗｐｓｕｂ］、［ｔｅｎｓｏｒＸｐｓｕｂ］の算出においても同様に全レコードをまとめて処理している。
　ステップ１２～１４において、演算部１２０は、各時点のΣｊ∈Ｒｉ→ｚｊｅｘｐ（→βＴ→ｚｊ）に対応する［Ｗｐｓｕｂ］を算出する。ステップ１６～１８において、演算部１２０は、各時点のΣｊ∈Ｒｉ→ｚｊ→ｚｊ Ｔｅｘｐ（→βＴ→ｚｊ）に対応する［ｔｅｎｓｏｒＸｐｓｕｂ］を算出する。
　ステップ２０において、演算部１２０は、Σｊ∈Ｒｉｅｘｐ（→βＴ→ｚｊ）の逆数に対応する［→ｙ］を算出する。逆数の計算はこの部分だけである。
　ステップ２２～２５において、演算部１２０は、勾配である式（５）を計算する。［→ｙ］、［→ｄ］は長さが時点数のベクトルであり、［Ｇ］、［Ｗ］、［Ｓ］は時点数×特徴量数の行列であり、計算結果の［ｔ→ｇ］は長さが特徴量数の横ベクトルである。ステップ２５のｓｕｍにより、全時点の総和を計算する。
　ステップ２７～３０において、演算部１２０は、へシアンである式（６）を計算する。ステップ２９の［Ｇｔｍｐ］［Ｇｔｍｐ］Ｔが、前述したＡＡＴに対応している。
　（実施の形態の効果）
　以上説明した本実施の形態に係る技術により、時点の数を復号することなく、効率良くｃｏｘ比例ハザードモデルのパラメータ推定を行うことが可能となる。
　すなわち、本実施の形態に係る技術により、従来の平文での処理のような反復処理を行わず全データを一度に処理することで、時点数を復号せずに計算できるようになる。また、反復処理を減らしたことにより、除算、指数、ｇｒｏｕｐ－ｂｙ　ｓｕｍと言った秘密計算での処理コストが大きい処理が並列化され、データや時点数を秘匿したまま効率良くｃｏｘ比例ハザード回帰のパラメータ推定が行えるようになる。
　（実施の形態のまとめ）
　本明細書には、少なくとも下記の各項に記載したパラメータ推定装置、パラメータ推定システム、パラメータ推定方法、及びプログラムが記載されている。
（第１項）
　ｃｏｘ比例ハザードモデルのパラメータ推定を、秘密計算により実行するパラメータ推定装置であって、
　イベントが観測された時点と、当該時点の観測対象の特徴量と、当該時点の観測対象の状態とを含むレコードを、観測対象毎に有するデータベースを格納するデータ格納部と、
　前記データベースから、時点からなるベクトルを読み出し、当該ベクトルをソートすることにより、置換表と、時点の境目を示すフラグとを生成し、前記置換表と、前記フラグとを用いることにより、時点の値を秘匿したまま、前記特徴量の時点毎の集計を行い、集計結果に基づいて前記パラメータ推定を行う演算部と、
　前記演算部により推定されたパラメータを出力する出力部と
　を備えるパラメータ推定装置。
（第２項）
　前記演算部は、前記パラメータ推定のための反復計算において使用される計算式における複数のｅｘｐの計算を、１反復あたり１回のｅｘｐの計算と、その計算結果を利用した演算により実行する
　第１項に記載のパラメータ推定装置。
（第３項）
　前記演算部は、前記パラメータ推定のための反復計算において使用される計算式における複数の逆数の計算を、１反復あたり１回の逆数の計算と、その計算結果を利用した演算により実行する
　第１項又は第２項に記載のパラメータ推定装置。
（第４項）
　前記演算部は、前記パラメータ推定のための反復計算において使用される計算式における時点毎の計算を、ベクトル又は行列又はテンソルを利用して全時点分についてまとめて実行する
　第１項ないし第３項のうちいずれか１項に記載のパラメータ推定装置。
（第５項）
　ｃｏｘ比例ハザードモデルのパラメータ推定を、秘密計算により実行するパラメータ推定システムであって、
　イベントが観測された時点と、当該時点の観測対象の特徴量と、当該時点の観測対象の状態とを含むレコードを、観測対象毎に有するデータベースを格納するデータ格納部と、
　前記データベースから、時点からなるベクトルを読み出し、当該ベクトルをソートすることにより、置換表と、時点の境目を示すフラグとを生成し、前記置換表と、前記フラグとを用いることにより、時点の値を秘匿したまま、前記特徴量の時点毎の集計を行い、集計結果に基づいて前記パラメータ推定を行う演算部と、
　前記演算部により推定されたパラメータを出力する出力部と
　を備えるパラメータ推定システム。
（第６項）
　ｃｏｘ比例ハザードモデルのパラメータ推定を、秘密計算により実行するパラメータ推定装置により実行されるパラメータ推定方法であって、
　イベントが観測された時点と、当該時点の観測対象の特徴量と、当該時点の観測対象の状態とを含むレコードを、観測対象毎に有するデータベースから、時点からなるベクトルを読み出し、当該ベクトルをソートすることにより、置換表と、時点の境目を示すフラグとを生成し、前記置換表と、前記フラグとを用いることにより、時点の値を秘匿したまま、前記特徴量の時点毎の集計を行い、集計結果に基づいて前記パラメータ推定を行う演算ステップと、
　前記演算ステップにより推定されたパラメータを出力する出力ステップと
　を備えるパラメータ推定方法。
（第７項）
　コンピュータを、第１項ないし第４項のうちいずれか１項に記載のパラメータ推定装置における各部として機能させるためのプログラム。
　以上、本実施の形態について説明したが、本発明はかかる特定の実施形態に限定されるものではなく、特許請求の範囲に記載された本発明の要旨の範囲内において、種々の変形・変更が可能である。
１００　パラメータ推定装置
１１０　入力部
１２０　演算部
１３０　出力部
１４０　格納部
１０００　ドライブ装置
１００１　記録媒体
１００２　補助記憶装置
１００３　メモリ装置
１００４　ＣＰＵ
１００５　インタフェース装置
１００６　表示装置
１００７　入力装置
１００８　出力装置

Claims (7)

  1. 　ｃｏｘ比例ハザードモデルのパラメータ推定を、秘密計算により実行するパラメータ推定装置であって、
    　イベントが観測された時点と、当該時点の観測対象の特徴量と、当該時点の観測対象の状態とを含むレコードを、観測対象毎に有するデータベースを格納するデータ格納部と、
    　前記データベースから、時点からなるベクトルを読み出し、当該ベクトルをソートすることにより、置換表と、時点の境目を示すフラグとを生成し、前記置換表と、前記フラグとを用いることにより、時点の値を秘匿したまま、前記特徴量の時点毎の集計を行い、集計結果に基づいて前記パラメータ推定を行う演算部と、
    　前記演算部により推定されたパラメータを出力する出力部と
    　を備えるパラメータ推定装置。
  2. 　前記演算部は、前記パラメータ推定のための反復計算において使用される計算式における複数のｅｘｐの計算を、１反復あたり１回のｅｘｐの計算と、その計算結果を利用した演算により実行する
    　請求項１に記載のパラメータ推定装置。
  3. 　前記演算部は、前記パラメータ推定のための反復計算において使用される計算式における複数の逆数の計算を、１反復あたり１回の逆数の計算と、その計算結果を利用した演算により実行する
    　請求項１又は２に記載のパラメータ推定装置。
  4. 　前記演算部は、前記パラメータ推定のための反復計算において使用される計算式における時点毎の計算を、ベクトル又は行列又はテンソルを利用して全時点分についてまとめて実行する
    　請求項１ないし３のうちいずれか１項に記載のパラメータ推定装置。
  5. 　ｃｏｘ比例ハザードモデルのパラメータ推定を、秘密計算により実行するパラメータ推定システムであって、
    　イベントが観測された時点と、当該時点の観測対象の特徴量と、当該時点の観測対象の状態とを含むレコードを、観測対象毎に有するデータベースを格納するデータ格納部と、
    　前記データベースから、時点からなるベクトルを読み出し、当該ベクトルをソートすることにより、置換表と、時点の境目を示すフラグとを生成し、前記置換表と、前記フラグとを用いることにより、時点の値を秘匿したまま、前記特徴量の時点毎の集計を行い、集計結果に基づいて前記パラメータ推定を行う演算部と、
    　前記演算部により推定されたパラメータを出力する出力部と
    　を備えるパラメータ推定システム。
  6. 　ｃｏｘ比例ハザードモデルのパラメータ推定を、秘密計算により実行するパラメータ推定装置により実行されるパラメータ推定方法であって、
    　イベントが観測された時点と、当該時点の観測対象の特徴量と、当該時点の観測対象の状態とを含むレコードを、観測対象毎に有するデータベースから、時点からなるベクトルを読み出し、当該ベクトルをソートすることにより、置換表と、時点の境目を示すフラグとを生成し、前記置換表と、前記フラグとを用いることにより、時点の値を秘匿したまま、前記特徴量の時点毎の集計を行い、集計結果に基づいて前記パラメータ推定を行う演算ステップと、
    　前記演算ステップにより推定されたパラメータを出力する出力ステップと
    　を備えるパラメータ推定方法。
  7. 　コンピュータを、請求項１ないし４のうちいずれか１項に記載のパラメータ推定装置における各部として機能させるためのプログラム。
PCT/JP2020/039119 2020-10-16 2020-10-16 パラメータ推定装置、パラメータ推定システム、パラメータ推定方法、及びプログラム WO2022079904A1 (ja)

Priority Applications (6)

Application Number Priority Date Filing Date Title
EP20957733.7A EP4231272A4 (en) 2020-10-16 2020-10-16 PARAMETER ESTIMATION DEVICE, PARAMETER ESTIMATION SYSTEM, PARAMETER ESTIMATION METHOD, AND PROGRAM
JP2022556813A JP7456514B2 (ja) 2020-10-16 2020-10-16 パラメータ推定装置、パラメータ推定システム、パラメータ推定方法、及びプログラム
US18/245,195 US20230367846A1 (en) 2020-10-16 2020-10-16 Parameter estimation apparatus, parameter estimation system, parameter estimation method, and program
AU2020472128A AU2020472128B2 (en) 2020-10-16 2020-10-16 Parameter estimation device, parameter estimation system, parameter estimation method, and program
CN202080106085.5A CN116324935A (zh) 2020-10-16 2020-10-16 参数估计装置、参数估计系统、参数估计方法及程序
PCT/JP2020/039119 WO2022079904A1 (ja) 2020-10-16 2020-10-16 パラメータ推定装置、パラメータ推定システム、パラメータ推定方法、及びプログラム

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
PCT/JP2020/039119 WO2022079904A1 (ja) 2020-10-16 2020-10-16 パラメータ推定装置、パラメータ推定システム、パラメータ推定方法、及びプログラム

Publications (1)

Publication Number Publication Date
WO2022079904A1 true WO2022079904A1 (ja) 2022-04-21

Family

ID=81209023

Family Applications (1)

Application Number Title Priority Date Filing Date
PCT/JP2020/039119 WO2022079904A1 (ja) 2020-10-16 2020-10-16 パラメータ推定装置、パラメータ推定システム、パラメータ推定方法、及びプログラム

Country Status (6)

Country Link
US (1) US20230367846A1 (ja)
EP (1) EP4231272A4 (ja)
JP (1) JP7456514B2 (ja)
CN (1) CN116324935A (ja)
AU (1) AU2020472128B2 (ja)
WO (1) WO2022079904A1 (ja)

Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
WO2006087854A1 (ja) * 2004-11-25 2006-08-24 Sharp Kabushiki Kaisha 情報分類装置、情報分類方法、情報分類プログラム、情報分類システム
WO2017119211A1 (ja) * 2016-01-07 2017-07-13 ソニー株式会社 情報処理装置、情報処理システム、および情報処理方法、並びにプログラム

Family Cites Families (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP6009697B2 (ja) * 2014-01-17 2016-10-19 日本電信電話株式会社 秘密計算方法、秘密計算システム、ソート装置及びプログラム

Patent Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
WO2006087854A1 (ja) * 2004-11-25 2006-08-24 Sharp Kabushiki Kaisha 情報分類装置、情報分類方法、情報分類プログラム、情報分類システム
WO2017119211A1 (ja) * 2016-01-07 2017-07-13 ソニー株式会社 情報処理装置、情報処理システム、および情報処理方法、並びにプログラム

Non-Patent Citations (1)

* Cited by examiner, † Cited by third party
Title
D.R.COX: "Regression Models and Life-Tables", JOURNAL OF THE ROYAL STATISTICAL SOCIETY. SERIES B (METHODOLOGICAL, vol. 34, no. 2, 1972, pages 187 - 220

Also Published As

Publication number Publication date
AU2020472128B2 (en) 2023-11-30
AU2020472128A1 (en) 2023-05-11
JPWO2022079904A1 (ja) 2022-04-21
CN116324935A (zh) 2023-06-23
EP4231272A1 (en) 2023-08-23
US20230367846A1 (en) 2023-11-16
EP4231272A4 (en) 2024-07-10
JP7456514B2 (ja) 2024-03-27

Similar Documents

Publication Publication Date Title
Blatt et al. Optimized homomorphic encryption solution for secure genome-wide association studies
Ablinger et al. The two-mass contribution to the three-loop gluonic operator matrix element Agg, Q (3)
Gu et al. A division-free Toom–Cook multiplication-based Montgomery modular multiplication
Borinsky et al. Tropical Feynman integration in the Minkowski regime
CN109328377B (zh) 秘密计算系统、秘密计算装置、秘密计算方法、以及程序
de Panafieu Phase transition of random non-uniform hypergraphs
CN112005288B (zh) 秘密聚合中值系统、秘密计算装置、秘密聚合中值方法、以及记录介质
WO2022079904A1 (ja) パラメータ推定装置、パラメータ推定システム、パラメータ推定方法、及びプログラム
JP6585846B2 (ja) 秘密計算システム、秘密計算装置、秘密計算方法、およびプログラム
Rossi Direct sampling of the self-energy with Connected Determinant Monte Carlo
JP7091930B2 (ja) テンソルデータ計算装置、テンソルデータ計算方法及びプログラム
US10333697B2 (en) Nondecreasing sequence determining device, method and program
Mounica et al. Implementation of 5-Qubit approach-based Shor's Algorithm in IBM Qiskit
Denisenko Quantum circuits for S-box implementation without ancilla qubits
Horacsek et al. A closed PP form of box splines via Green’s function decomposition
JP2022083992A (ja) 暗号化ハッシュ関数の原像要素を決定するための方法、コンピュータプログラム、およびデータ処理システム
Takahashi JMASM44: Implementing Multiple Ratio Imputation by the EMB Algorithm (R)
Xu et al. Efficient algorithms for computing multidimensional integral fractional Laplacians via spherical means
Ruffa et al. Parallelized solution of banded linear systems with an introduction to p-adic computation
WO2022254599A1 (ja) 秘密共役勾配法計算方法、秘密共役勾配法計算システム、秘密計算装置、およびプログラム
JP7173328B2 (ja) 秘密除算システム、秘密計算装置、秘密除算方法、およびプログラム
Petković et al. Ostrowski-like method with corrections for the inclusion of polynomial zeros
Miller Quantum resource counts for operations constructed from an addition circuit
WO2024079899A1 (ja) 統計値推定装置、統計値推定システム、統計値推定方法、及びプログラム
WO2022124010A1 (ja) 演算制御装置、演算制御方法、および記録媒体

Legal Events

Date Code Title Description
121 Ep: the epo has been informed by wipo that ep was designated in this application

Ref document number: 20957733

Country of ref document: EP

Kind code of ref document: A1

ENP Entry into the national phase

Ref document number: 2022556813

Country of ref document: JP

Kind code of ref document: A

WWE Wipo information: entry into national phase

Ref document number: 2020472128

Country of ref document: AU

ENP Entry into the national phase

Ref document number: 2020472128

Country of ref document: AU

Date of ref document: 20201016

Kind code of ref document: A

NENP Non-entry into the national phase

Ref country code: DE

ENP Entry into the national phase

Ref document number: 2020957733

Country of ref document: EP

Effective date: 20230516