Patent Application
20230367846
The present invention relates to a technique for implementing parameter estimation of a cox proportional hazard model while keeping confidentiality in a secure computation technique.
A cox proportional hazard regression which is a regression analysis using a proportional hazard model of cox is an analysis method (NPL 1) frequently used in survival time analysis. In commercial statistical analysis software, or a package of software language, cox proportional hazard regression can be performed in plain text.
A method called secure computation is known as a method to obtain a specific calculation result without restoring the encrypted numerical value. As an example, by encrypting numerical fragments by distributing them among several secure computation devices and having multiple secure computation devices perform cooperative computation, the results of addition, subtraction, constant addition, results of addition and subtraction, constant addition, multiplication, constant multiplication, logical operations (negation, logical product, logical sum, and exclusive disjunction), and data format conversion (integers and binary numbers), etc. can be obtained in a distributed manner in several secure computation devices, without restoring numerical values.
When the parameter estimation of the cox proportional hazard model is performed, calculation is repeated every point of time when death or truncation occurs in a plain text. However, in the case of using the same method in the secure computation, there is a problem that it is necessary to decode the value at the point of time that should be kept secret. In addition, since the calculation of parameter estimation of the cox proportional hazard model includes many processes such as an exponent, division, and group-by sum, which are high in processing cost, and it is difficult to efficiently calculate by secure computation.
The present invention has been made in view of the above-mentioned points, and an object of the present invention is to provide a technique for efficiently performing parameter estimation of a cox proportional hazard model without decoding values at the point of time.
According to the disclosed technique, a parameter estimation device, which executes a parameter estimation of a cox proportional hazard model by secure computation comprises:
According to the disclosed technique, there is provided a technique for efficiently performing parameter estimation of a cox proportional hazard model without decoding a value at a point of time.
Hereinafter, an embodiment of the present invention (the present embodiment) will be described with reference to the drawings. The embodiment to be described below is merely exemplary and embodiments to which the present invention is applied are not limited to the following embodiment.
The parameter estimation device 100 may be constituted by one device (computer), or constituted by a system consisting of a plurality of computers. This system may be called a parameter estimation system. In the parameter estimation system, for example, the calculation unit 120 and the data storage unit 140 may be separate servers.
Concealed data from data obtained by observation is input to an input unit 110 of the parameter estimation device 100. Unless otherwise stated, the data handled by the parameter estimation device 100 is concealed data, and calculation shall be performed by secure computation.
The input data is stored in the data storage unit 140 as a database. The calculation unit 120 performs parameter estimation of the cox proportional hazard model by performing processing to be described later on data such as scalar, vector, matrix and the like read from the database of the data storage unit 140. The output unit 130 outputs a parameter 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 in response to an external access. The processing details by 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 details of processing described in the present embodiment. Note that this “computer” may be a physical machine, or may be a virtual machine in the cloud. When a virtual machine is used, “hardware” described here is virtual hardware.
The program can be recorded on a computer-readable recording medium (a portable memory or the like) to store or distribute the program. It is also possible to provide the program through a network such as the Internet or e-mail.
A program for realizing processing in the computer is provided by, for example, a recording medium 1001 such as a CD-ROM or a memory card. When the recording medium 1001 having the program stored therein is set in the drive device 1000, the program is installed in the auxiliary storage device 1002 from the recording medium 1001 via the drive device 1000. However, the program does not necessarily have to be installed from the recording medium 1001, and may be downloaded from another computer via a 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 start instruction of the program is received. The CPU 1004 realizes a function related to the parameter estimation device 100 according to a 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 GUI (Graphical User Interface) and the like by the program. The input device 1007 is configured 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 a calculation result.
Before describing an example of operation of the parameter estimation device 100, notation, cox proportional hazard regression, secure computation and the like will be described as preparation.
Defining a by b is described as a:=b, vector is described as →a:=(a0, . . . , an−1), unless otherwise specified, matrix is described as capital letter such as A, and a transposed matrix of that matrix is described as AT. In the text of the present specification, for convenience of description, the vector symbol (→) to be placed above the beginning of the character is described before the character, such as “→a”. Further, in the algorithm described in the drawing, tensor of third floor is described in italic characters. In the text of the present specification, for convenience of description, characters of tensor of third floor are represented, for example, by attaching a tensor to the upper left of the characters, as shown by tensorZZ′.
In addition, subtraction and multiplication, when the input is a vector →a or a matrix A and a scalar b, it is assumed that the operation with b is performed on all the elements of →a and A. The vector not specifically described is a column vector. In the case of a row vector, it is distinguished by adding t at the upper left like t→a.
The cox proportional hazard model is a model represented by the equation (1) (NPL 1).
[Math. 1]
λ(t|{right arrow over (z)})=λ0(t)exp({right arrow over (β)}T{right arrow over (z)}) (1)
In equation (1), t, β, and z represent time, weight, and feature amount, respectively, and λ0(t) and exp(→βT→z) are called baseline hazard function and relative risk function (hazard), respectively. In cox proportional hazards regression, the weight parameter β is estimated. The weight is calculated by performing maximum likelihood estimation using the part likelihood shown in equation (2).
In equation (2), D is the number of a point of time when death was observed, and →zi represents the feature amount of the patient who died at point of time i. In this embodiment, death is used as an objective event, but this is an example. For example, the event may be for the purpose of falling, inability to walk, onset of illness, hospitalization, or the like. Truncation may be interpreted as one of the events.
Ri in equation (2) is a set of patients who have not been truncated or died until just before a point of time i, and is called a risk set. In addition, “truncation” means that observation becomes impossible and it becomes unknown whether or not death has occurred after that.
Therefore, the part likelihood function of equation (2) is calculated by calculating (hazard of deceased patient)/(sum of hazards of risk set) for each point of time and multiplying by all points of time. Since this part likelihood is based on the assumption that multiple truncations and deaths do not occur at the same time (no tie data), the actual data that often has tie data is the Breslow method shown in equation (3) is often used.
Equation (3) is basically the same as the part likelihood of cox in equation (2), but the denominator is di-th power (di is the number of dead patients at point of time i), and the difference is that the sum si of the feature amount of the dead patients at point of time i is used instead of →zi in the numerator. Since the case of di=1 agrees with the equation (2), the equation (3) can be considered as a generalized form of the equation (2). In the following description, equation (3) is premised.
The Newton method or the like is generally used as a method for obtaining the maximum likelihood estimation of →β. In this embodiment, the Newton method is used. In 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 (hessian) of the log-likelihood function. The log-likelihood function l (of →β), its first derivative U(of →β), and its second derivative I(of →β) are shown in equation (4), equation (5), and equation (6), respectively.
In Newton's method, the following equation (7) is repeated using equations (5) and (6) to obtain the maximum likelihood estimation value of →β. Equation (7) converges after about 5 iterations.
[Math. 7]
{right arrow over (β)}={right arrow over (β)}+I({right arrow over (β)})−1U({right arrow over (β)}) (7)
A value that is concealed a value of a by encryption, secret sharing or the like is called a ciphertext or concealed value, is described as [a]. When a is concealed by secret sharing, a set of secret sharing fragments, which is had by each secure computation device, is referred by the [a]. Further, the parenthesis “[“, ”]” 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 “[“, ”]”.
Next, various calculations in the secure computation will be described.
Addition, subtraction, and multiplication of two ciphertexts [a] and [b] are processes for calculating ciphertexts [a +b], [a−b], and [a×b], respectively. These calculations are described as [a]+[b], [a]−[b], and [a]×[b], respectively.
Further, the process of dividing the ciphertext [a] by the plaintext b uses a notation such as [a]/b. When the input is a vector or a matrix and these processes are performed for each element, the notation is similarly [→a]+[→b], [A]+[B].
When the input is the matrix a and the column vector →b, in the addition, subtraction and multiplication, the arithmetic operation of each element →b is performed for each column vector of the matrix, and when the input is the matrix a and the row vector t→b, the arithmetic operation of each element t→b is performed for each row vector of the matrix.
The process of finding the sum of the elements of the vector [→a] is described as sum ([→a]). When the matrix of m×n is sum input as in sum ([A]), the sum in the column direction is calculated and the row vector [t→c] of length n is output.
The calculation of the vector such as ([b], [b]−[a1], [b]−([a1]+[a2]), . . . , [b]−Σ[→a]) is described as prefixSub ([→a], [b]), from vector [→a]:=([a1], [a2], . . . , [an]) and scalar [b].
The calculation of the reciprocal 1/[a] of the ciphertext [a] is described as [c]←reciprocal([a]). The same notation is used when the input is a vector.
Entering the ciphertext [a] and calculating the [a] power of the Napier number e is described as [c]←exp([a]). The same notation is used when the input is a vector.
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 replacement table [→π] and a flag [→e] indicating whether the value of the key is a boundary or not, and by reusing these, various Group-by operations using the same key can be performed efficiently.
With reference to
The replacement table →π is a vector indicating how each element in the vector inputted to Group-by common move to what number of elements the element of the vector is sorted.
In
The flag →e is a vector in which 0 is input to the position of the element for each element of the sorted vector of the vector input to the Group-by common if the value is the same the value below the element by comparison, and 1 is input to the position of the element if the value is different from the value below the element.
The Group-by Common is performed by inputting the vector of the key [→k] is described as shown by the equation (8). The process of sorting the vector [→a] and the matrix [A] (the number of rows is equal to the length of [→π]), using the replacement table [→π]is described as in equations (9) and (10).
It is described as in equations (11) and (12) that Group-by sum is performed using the sorted vectors [→a′] and [A′] and [→e]. Primes such as [→a′] and [A′] indicate that they have been sorted, and the same notation is used thereafter.
[Math. 8]
{right arrow over (π)}
,
{right arrow over (e)}
←groupByCommon(
{right arrow over (k)}
) (8)
[Math. 9]
{right arrow over (a′)}
←sort(
{right arrow over (a)}
,
{right arrow over (π)}
) (9)
[Math. 10]
A′←sort(
A
,
{right arrow over (π)}
) (10)
[Math. 11]
{right arrow over (c)}
←groupBySum(
{right arrow over (a′)}
,
{right arrow over (e)}
) (11)
[Math. 12]
C
←groupBySum(
A′
,
{right arrow over (e)}
) (12)
In the case of the input of sort or groupBySum is a matrix, the processing is performed for each column. Although the size of the output is generally equal to or less than the size of the input by performing Group-by sum, in the present embodiment, the sizes of the input and the output are the same, and the unnecessary portion is padded at the tail with zero. This makes it possible to conceal the number of attributes of key of point of time. In the explanation of the process described later, for convenience, the result of Group-by sum is described as “a vector whose length is the number of points of time” or “a matrix of the number of points of time×n”, in fact, to correspond to “a vector whose length is the number of records” and “a matrix with the number of records×n”, the unnecessary part is padded to 0.
For example, assuming that the flag →e shown in
The operation example of the parameter estimation device 100 will be described below. 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 equation (5), equation (6), and equation (7) by secure computation, and then, the parameters of the cox proportional hazard regression are estimated. In the following, first, the characteristic operation will be described.
When the equations (5) and (6) are mounted as they are, the calculation is repeated at each point of time, and the calculation results are added in order. In the case of calculating the secure computation cox proportional hazard regression, previously mentioned groupByCommon and groupBySum are used in order to keep the number of points of time secret.
In groupByCommon and groupBySum, totalization is performed for each point of time while keeping the value of the point of time secret, and since unnecessary parts are padded at zero, information on the number of points of time is not leaked.
Further, by processing all the records collectively instead of repeating the processing for each number of points of time, the number of operations having a large processing cost in the secure computation can be reduced and the processing efficiency is good. That is, Σj∈Ri exp(→βT→zj) for all points of time is calculated at once, instead of calculation, for example, that Σj∈Ri exp(→βT→zj) is calculated at point of time 1, and Σj∈Ri exp(→βT→zj) is calculated at point of time 2, and . . . . More specifically, when the value at each point of time is a scalar, the length is treated as a vector of the number of points of time and calculated collectively, and when the value at each point of time is a vector with length n, it is treated and calculated all together as a matrix of the number of points of time×n. If the value at each point of time is an n×n matrix, it is treated as a tensor of third floor with the number of points of time×n×n and calculated collectively.
Equation (5) and equation (6) include many exponentiations and divisions, and since the processing of Σj∈Ri is Group-by sum, the calculation cost in secure computation is very large.
The parameter estimation device 100 according to the present embodiment minimizes costly processing such as exp, division, and Group-by sum in cox proportional hazard regression, and calculates efficiently. When simply calculated according to equation (5) and equation (6), exponentiation is required for 7 times×number of points of time and division is required for 3 times×number of points of time per iteration of Newton's method, but in this embodiment, the number of points of time was minimized as follows.
As shown in the expression (5) and expression (6), since all the arguments of EXP are →βT→ZJ, it is sufficient to reuse the arguments once calculated. In addition to this, processing is made parallel by collectively processing all the records, and it is possible to complete the processing once per iteration.
When division is performed by reciprocal calculation multi plication, the reciprocal of Σj∈Ri exp(→βT→zj) can be reused in the second term of equation (5) and the first term of equation (6), and only need to calculate the once reciprocal+twice multiplication instead of dividing twice. The second term of equation (6) has a different divisor from the above two, but this term can be obtained without division. Comparing the second term of the formula (5) and the second term of the formula (6), when the part excluding the di of the second term of the formula (5) is A, the second term of the formula (6) Since the term can be expressed by AAT, the second term of equation (6) can be calculated only by the product. In addition to this, by processing all the above-mentioned records at once, the processing is parallelized, and only once per iteration is required.
If implemented according to equation (5) and equation (6), Group-by sum will be performed many times during processing, and processing efficiency will decrease. Therefore, in the present embodiment, paying attention to the fact that all the keys are the same, it is decided to reuse [→e] obtained by Group-by common which was performed only once for the key at the beginning.
In the parameter estimation device 100 of the present embodiment, the Group-by common is utilized to efficiently calculate the secure computation cox proportional hazard regression. In the calculation of groupBySum, only the totalization using the flag [→e] indicating the boundary is performed.
Next, details of the processing executed by the parameter estimation device 100 will be described. Here, it is assumed that the concealed observation data is stored as a database in a data storage unit 140, and a calculation unit 120 estimates a parameter of the cox proportional hazard regression by processing the data by secure computation. In the processing operation, the above-mentioned characteristic operation is performed.
As shown in
In the example of
At the same time, death or truncation may be observed for a plurality of patients. Therefore, the number of points of time D≤the number of patients m.
A calculation unit 120 reads the data from the database to hold the feature amount of all the patients as a matrix Z of M×N, a point of time as a time vector →t, and a state of all the patients as a state vector →C.
None of Z, →t, and →c is sorted in the initial state. As mentioned above, by Group-by common, the time vector →t is first sorted to create a replacement table →π, and by reusing it, and the sorting by the key of the feature amount Z or →t of state vector →c requires only sorting 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 estimates parameters of the data in the database according to the procedure of the algorithm shown in
In step 3 of the algorithm 1 of
In steps 5 and 6, the calculation unit 120 sorts [Z] and [→c], respectively, to create [Z′] and [→c′].
in step 8, the calculation unit 120 creates [Z′dead] which sets all feature amounts other than the dead cases feature amounts from [Z′] to 0, and in step 9, [S] of the sum of the features of the dead cases at each point of time by groupBySum is created. In step 11, the number of deaths [→d] at each point of time is calculated. In step 13, [tensorZZ′], which is an m×n×n tensor, is created, and [→β] is updated in steps 15 to 17.
The processing of calcGH in step 16 will be described with reference to
In steps 8 to 10, the calculation unit 120 calculates [→vpsub] corresponding to Σj∈Ri exp(→βT→zj) at each point of time. Each element of [→vpsub] , which is a vector of length of the number of points of time, has Σj∈Ri exp(→βT→zj) in a scalar value. That is, the calculation here is not an iterative process for each number of points of time, but a calculation that processes all records at once. In the calculation of [Wpsub] and [tensorXpsub] described below, all the records are processed together in the same manner.
In steps 12 to 14, the calculation unit 120 calculates [Wpsub] corresponding to Σj∈Ri→zj exp(→βT→zj) at each point of time. In steps 16 to 18, the calculation unit 120 calculates [tensorXpsub] corresponding to j∈Ri→zj exp(→βT→zj) at each point of time.
In step 20, the calculation unit 120 calculates [→y] corresponding to the reciprocal of Σj∈Ri exp(→βT→zj). This is the only part of the reciprocal calculation.
In steps 22 to 25, the calculation unit 120 calculates the equation (5) which is a gradient. [→y] and [→d] are vectors whose length is the number of points of time, [G], [W] and [S] are matrices of the number of points of time multiplexed by the number of feature amounts, and the calculation result [t→g] is a column vector whose length is the number of features. The sum of all points of time is calculated by sum in step 25.
In steps 27 to 30, the calculation unit 120 calculates the equation (6) which is a hessian. [Gtmp][Gtmp]T in step 29 corresponds to the above-mentioned AAT.
According to the technique according to the present embodiment described above, it is possible to efficiently perform parameter estimation of the cox proportional hazard model without decoding the number of points of time.
That is, by processing all data at a time without performing repetitive processing such as processing in a conventional plaintext according to the technique of the present embodiment, the number of points of time can be calculated without decoding. Further, by reducing the repetition processing, processing such as division, exponent, and group-by sum which is high in processing cost in secure computation is parallelized, and the parameter estimation of the cox proportional hazard regression can be efficiently performed while keeping the data and the number of points of time secret.
At least the parameter estimation device, the parameter estimation system, the parameter estimation method, and the program are described in the present specification.
A parameter estimation device, which executes a parameter estimation of a cox proportional hazard model by secure computation, a parameter estimation device comprising:
The parameter estimation device according to item 1, wherein the calculation unit executes a calculation of a plurality of exps in a calculation equation used in an iterative calculation for the parameter estimation by the calculation of exp once per iteration and a calculation using the calculation result.
The parameter estimation device according to item 1 or 2, wherein
The parameter estimation device according to one of items 1 to 3, wherein
A parameter estimation system, which executes a parameter estimation of a cox proportional hazard model by secure computation, a parameter estimation system comprising:
A parameter estimation method executed by a parameter estimation device that executes a parameter estimation of a cox proportional hazard model by secure computation, a parameter estimation method comprising:
The program that causes a computer to function as each unit of the parameter estimation device according to any one of items 1 to 4.
The embodiment has been described above, but the present invention is not limited to the specific embodiment. Various modifications and changes can be made within the scope of the gist of the present invention described in the claims.
| Filing Document | Filing Date | Country | Kind |
|---|---|---|---|
| PCT/JP2020/039119 | 10/16/2020 | WO |
