WO2018168383A1 - 最適解判定方法、最適解判定プログラム及び最適解判定装置 - Google Patents

最適解判定方法、最適解判定プログラム及び最適解判定装置 Download PDF

Info

Publication number
WO2018168383A1
WO2018168383A1 PCT/JP2018/006501 JP2018006501W WO2018168383A1 WO 2018168383 A1 WO2018168383 A1 WO 2018168383A1 JP 2018006501 W JP2018006501 W JP 2018006501W WO 2018168383 A1 WO2018168383 A1 WO 2018168383A1
Authority
WO
WIPO (PCT)
Prior art keywords
solution
evaluation value
solutions
maximum
optimal solution
Prior art date
Application number
PCT/JP2018/006501
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 JP2019505816A priority Critical patent/JP6851460B2/ja
Priority to CN201880009956.4A priority patent/CN110249344A/zh
Priority to EP18768409.7A priority patent/EP3598350A4/en
Publication of WO2018168383A1 publication Critical patent/WO2018168383A1/ja
Priority to US16/513,202 priority patent/US11816580B2/en

Links

Images

Classifications

    • G—PHYSICS
    • G16—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR SPECIFIC APPLICATION FIELDS
    • G16B—BIOINFORMATICS, i.e. INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR GENETIC OR PROTEIN-RELATED DATA PROCESSING IN COMPUTATIONAL MOLECULAR BIOLOGY
    • G16B40/00—ICT specially adapted for biostatistics; ICT specially adapted for bioinformatics-related machine learning or data mining, e.g. knowledge discovery or pattern finding
    • 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/11—Complex mathematical operations for solving equations, e.g. nonlinear equations, general mathematical optimization problems
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06N—COMPUTING ARRANGEMENTS BASED ON SPECIFIC COMPUTATIONAL MODELS
    • G06N5/00—Computing arrangements using knowledge-based models
    • G06N5/01—Dynamic search techniques; Heuristics; Dynamic trees; Branch-and-bound
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06Q—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR ADMINISTRATIVE, COMMERCIAL, FINANCIAL, MANAGERIAL OR SUPERVISORY PURPOSES; SYSTEMS OR METHODS SPECIALLY ADAPTED FOR ADMINISTRATIVE, COMMERCIAL, FINANCIAL, MANAGERIAL OR SUPERVISORY PURPOSES, NOT OTHERWISE PROVIDED FOR
    • G06Q10/00—Administration; Management
    • G06Q10/04—Forecasting or optimisation specially adapted for administrative or management purposes, e.g. linear programming or "cutting stock problem"
    • G—PHYSICS
    • G16—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR SPECIFIC APPLICATION FIELDS
    • G16B—BIOINFORMATICS, i.e. INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR GENETIC OR PROTEIN-RELATED DATA PROCESSING IN COMPUTATIONAL MOLECULAR BIOLOGY
    • G16B25/00—ICT specially adapted for hybridisation; ICT specially adapted for gene or protein expression
    • G16B25/10—Gene or protein expression profiling; Expression-ratio estimation or normalisation
    • G—PHYSICS
    • G16—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR SPECIFIC APPLICATION FIELDS
    • G16B—BIOINFORMATICS, i.e. INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR GENETIC OR PROTEIN-RELATED DATA PROCESSING IN COMPUTATIONAL MOLECULAR BIOLOGY
    • G16B30/00—ICT specially adapted for sequence analysis involving nucleotides or amino acids
    • G—PHYSICS
    • G16—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR SPECIFIC APPLICATION FIELDS
    • G16B—BIOINFORMATICS, i.e. INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR GENETIC OR PROTEIN-RELATED DATA PROCESSING IN COMPUTATIONAL MOLECULAR BIOLOGY
    • G16B5/00—ICT specially adapted for modelling or simulations in systems biology, e.g. gene-regulatory networks, protein interaction networks or metabolic networks
    • G16B5/20—Probabilistic models
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06N—COMPUTING ARRANGEMENTS BASED ON SPECIFIC COMPUTATIONAL MODELS
    • G06N3/00—Computing arrangements based on biological models
    • G06N3/12—Computing arrangements based on biological models using genetic models
    • G06N3/126—Evolutionary algorithms, e.g. genetic algorithms or genetic programming

Definitions

  • the present invention relates to an optimal solution determination method, an optimal solution determination program, and an optimal solution determination device, and more particularly to a technique for determining the optimality of a solution in a combinatorial optimization problem.
  • the combinatorial optimization problem includes many difficult problems that are NP (Non-deterministic Polynomial time) complete or difficult to NP. That is, generally, as the problem scale increases, the amount of calculation explodes in the order of an exponential function or higher, so that a solution by exhaustive full search is almost impossible.
  • NP Non-deterministic Polynomial time
  • Non-Patent Document 1 In order to respond directly to this problem, there is a study to determine the optimality of the solution by assigning some score to the solution and using statistical analysis together, for example, the method of estimating the maximum score by extreme value statistical analysis, etc. (Non-Patent Document 1).
  • ZDD Zero-suppressed Binary Decision Diagram
  • RNA ribonucleic acid
  • NGS Next Generation Sequencer
  • Patent Document 1 discloses a technique using a Bayesian model based on nonparametric regression for estimating a network relationship between genes from a time series study of gene expression.
  • the drug discovery project is intended to produce drugs that can be used by the human body, and it takes a long time and enormous costs to clarify the results. high. Therefore, if the optimality of the graph structure that skillfully explains the acquired RNA expression matrix data can be determined, it is very useful in industry.
  • the method for constructing a gene regulatory network described in Patent Document 1 obtains time-dependent data on gene expression for several genes, determines the correction of Bayesian estimation methods, and the relationship between causes and results among expressed genes. To do. Modification of the Bayesian estimation method includes estimating causal relationships between expressed genes using time-varying data. In the method described in Patent Document 1, Bayesian estimation and nonparametric regression are corrected based on the assumption that a change in a slow-expressing gene is not likely to cause a change in a fast-expressing gene. While providing a high network solution, there is no guarantee that the above assumptions can be applied to all gene regulatory networks.
  • the present invention has been made in view of such circumstances, and an optimal solution determination method, an optimal solution determination program, and an optimal solution determination device capable of efficiently and accurately determining the optimality of a solution in a combinatorial optimization problem
  • the purpose is to provide.
  • an invention is an optimal solution determination method for determining optimality of a solution in a combination optimization problem by a computer, and determining the optimality of a solution in a combination optimization problem by a computer.
  • a first step of uniformly extracting a plurality of solutions on the solution space of the combinatorial optimization problem as a first plurality of solutions, and a first step of uniformly extracting by the first step A second step of acquiring a first plurality of evaluation values corresponding to each of the plurality of solutions, and a number exceeding the number of the first plurality of solutions based on the acquired first plurality of evaluation values
  • a plurality of solutions on the solution space of the combinatorial optimization problem are uniformly extracted as the first plurality of solutions, and the number exceeds the number of the first plurality of uniformly extracted solutions. Is estimated as the first maximum evaluation value. Then, when at least one of the solutions belonging to the solution space is acquired as a solution candidate, an evaluation value corresponding to the solution candidate is acquired, and whether the acquired evaluation value falls within the confidence interval of the first maximum evaluation value It is determined whether or not (whether the solution candidate is one of the first optimum solutions). Thereby, it is possible to determine the statistical optimality of the solution of the combination optimization problem, and it is possible to determine the optimality of whether or not the solution satisfies a sufficient condition.
  • the first plurality of solutions are U ⁇ V solutions, where U and V are natural numbers, respectively, and the third step is U ⁇ V solutions.
  • the V solution is divided into V blocks, and the V maximum value of the U solution evaluation values is obtained for each block, and the maximum value of the segment follows the general extreme value distribution using the V maximum values. It is preferable to estimate the first maximum evaluation value. Get the maximum number of partition values of evaluation values of U solutions for each block, and use the maximum number of partition values of V to determine that the maximum partition value follows a generalized extreme value distribution (GEV).
  • the evaluation value (first maximum evaluation value) is estimated by maximum likelihood.
  • the method further includes an eighth step of outputting a determination result in the seventh step.
  • a first search method having a low calculation cost but a low solution accuracy and a first calculation method having a higher calculation cost but a higher solution accuracy than the first search method.
  • the first step is to input a first solution candidate first searched by the first search method, and the evaluation value of the first solution candidate is the first maximum evaluation value. It is preferable to input the second solution candidate searched by the second search method only when it does not fall within the confidence interval.
  • the sufficiency determination fails. Since the heuristic search by the first search method is insufficient, it is suggested that another optimal solution exists. In this case, the search is switched to the search for the second solution candidate by the second search method with high calculation cost but high solution accuracy, and the sufficiency of the second solution candidate searched by the second search method is determined.
  • a plurality of solutions that are solutions in the solution space and whose distance in the solution space from the first optimum solution is out of a certain range are determined as the second plurality of solutions.
  • the tenth step of obtaining the second plurality of evaluation values corresponding to each of the second plurality of solutions, and the second plurality of evaluation values The maximum evaluation value when the number of solutions exceeding the number of the second plurality of solutions is assumed is estimated, and the eleventh step to be the second maximum evaluation value and the evaluation value of the first optimal solution are acquired. It is preferable to include a twelfth step and a thirteenth step of determining whether or not the evaluation value of the first optimal solution exceeds the second maximum evaluation value.
  • a plurality of solutions whose distance in the solution space from the first optimal solution is outside a certain range are uniformly extracted as the second plurality of solutions, and the second extracted uniformly.
  • the maximum evaluation value when the number of solutions exceeding the number of the plurality of solutions is assumed is estimated as the second maximum evaluation value.
  • the optimal solution of 1 satisfies the necessary and sufficient conditions and is the optimal solution in the solution space, and there is no other equivalent solution.
  • the second plurality of solutions are P ⁇ Q solutions where P and Q are natural numbers, respectively, and the eleventh step is P ⁇ Q.
  • the solution is divided into Q blocks, and the maximum number of pieces of evaluation values of the P pieces of solutions for each block is obtained.
  • the maximum number of pieces follows the general extreme value distribution using the maximum number of pieces of Q pieces. It is preferable to estimate the second maximum evaluation value. Get the maximum number of partition values of evaluation values of P solutions for each block, and use the maximum number of partition values of Q to determine that the maximum value of the partition follows the general extreme value distribution (second maximum evaluation) Value).
  • the method further includes a fourteenth step of outputting a determination result in the thirteenth step.
  • the certain range is expanded.
  • a fifteenth step of uniformly extracting a plurality of solutions outside the expanded fixed range as a third plurality of solutions and a third plurality of evaluation values corresponding to each of the third plurality of solutions Based on the 16 steps and the third plurality of evaluation values, the maximum evaluation value when the number of solutions exceeding the number of the third plurality of solutions is assumed is estimated to be the third maximum evaluation value. It is preferable to further include 17 steps and an 18th step of determining whether or not the evaluation value of the first optimal solution exceeds the third maximum evaluation value.
  • the evaluation value of the first optimal solution searched as one of the optimal solutions does not exceed the second maximum evaluation value, an equivalent optimal solution may exist.
  • the first optimal solution may exist.
  • the third maximum evaluation value is estimated again based on the third plurality of solutions uniformly extracted from the solution space away from the optimal solution (the solution space outside the expanded fixed range). Thereby, it can be confirmed that there is no equivalent optimum solution in the solution space away from the searched optimum solution.
  • the solution in the solution space is determined.
  • a nineteenth step of obtaining a fourth solution candidate that is a fixed distance away from the first optimum solution, and solutions from each of the first optimum solution and the fourth solution candidate in the solution space A twentieth step for uniformly extracting a fourth plurality of solutions whose spatial distances are outside a certain range, and a twenty-first evaluation for obtaining a fourth plurality of evaluation values corresponding to each of the fourth plurality of solutions.
  • step of preferably further comprises a.
  • the evaluation value of the first optimal solution searched as one of the optimal solutions does not exceed the second maximum evaluation value, the first optimal solution and a fourth distance that is a fixed distance away from the first optimal solution
  • a fourth plurality of solutions each having a distance in the solution space outside a certain range are uniformly extracted from both of the solution candidates, and a fourth maximum evaluation value is obtained based on the uniformly extracted fourth plurality of solutions. Estimate again. As a result, it can be confirmed that there is no equivalent optimal solution in the solution space apart from the searched first optimal solution and the fourth solution candidate outside the fixed range.
  • the solution space includes a first solution space in the first constraint condition and a second solution space in the second constraint condition, and the first solution
  • the maximum evaluation value in this case is estimated to be the fifth maximum evaluation value in the 27th step and the solution candidates in the second solution space, from the first solution candidate obtained in the 24th step
  • the solution when determining the optimality of a solution candidate in the first solution space under the first constraint condition, the solution is a nearby solution that is close in distance from the solution candidate in the solution space.
  • the neighborhood solution in the second solution space under the second constraint condition is acquired, and the optimality of the neighborhood solution is determined. That is, substituting the optimality determination of the neighborhood solution, and if the neighborhood solution is optimum, it is estimated that the solution candidate in the first solution space at a distance close to the neighborhood solution is also optimum.
  • the first step is inappropriate from the combinable patterns in the combinatorial optimization problem without considering the remaining combinations due to a part of the combinations.
  • the step of reducing the pattern to be identified By identifying whether or not is fixed, the step of reducing the pattern to be identified, and extracting the common part of the pattern group having a difference only in a part of the combination from the patterns that can be combined, and the remaining combinations Using a data structure that reduces and enumerates combinable patterns using at least one of the steps of reducing patterns to be identified by sharing, and calculates the total number of solutions in the solution space; It is preferable to generate a random number equal to or less than the calculated total number and extract a solution corresponding to a pattern specified by the generated random number.
  • the combination optimization problem is preferably a combination optimization problem of a gene regulatory network.
  • An optimal solution determination program causes a computer to execute the above-described optimal solution determination method.
  • the invention according to still another aspect is an optimal solution determination apparatus that determines the optimality of a solution in a combinational optimization problem, wherein a plurality of solutions on the solution space of the combinational optimization problem are defined as a first plurality of solutions.
  • the first evaluation value acquisition unit that acquires the evaluation value corresponding to each of the first plurality of uniformly extracted solutions, and the acquired first plurality of evaluation values,
  • a first maximum evaluation value estimating unit configured to estimate a maximum evaluation value when the number of solutions exceeding the number of the first plurality of solutions is assumed, and a first maximum evaluation value;
  • a solution acquisition unit that acquires at least one solution as a solution candidate, an evaluation value corresponding to the solution candidate is acquired from the first evaluation value acquisition unit, and the evaluation value corresponding to the acquired solution candidate is the first maximum evaluation value And whether the evaluation value corresponding to the solution candidate is the first maximum rating. If it is determined to fall within the confidence interval of values, comprising a first determination unit that the solution candidate to the first optimal solution, the.
  • a plurality of solutions that are solutions in the solution space and whose distance in the solution space from the first optimum solution is out of a certain range are determined as the second plurality of solutions.
  • a second solution extraction unit that uniformly extracts a solution
  • a second evaluation value acquisition unit that acquires a second plurality of evaluation values corresponding to each of the second plurality of solutions, and a second plurality of evaluations
  • a second maximum evaluation value estimating unit configured to estimate a maximum evaluation value when the number of solutions exceeding the number of the second plurality of solutions is assumed based on the value, and to set the second maximum evaluation value;
  • the evaluation value corresponding to the optimal solution is acquired from the second evaluation value acquisition unit, and the evaluation value corresponding to the acquired first optimal solution is determined whether or not it exceeds the second maximum evaluation value. It is preferable to further include two determination units.
  • Diagram showing RNA expression matrix data and diagram showing gene regulatory network Conceptual diagram showing features of the present invention
  • Functional block diagram showing a first embodiment of the present invention A diagram showing an entire set “G ⁇ ” composed of elements “A, B, C, D”.
  • FIGS. 1-10 A diagram showing four subsets (G (1) , G (2) , G (3) , G (4) ) The figure which shows the "pattern” corresponding to selection of three subsets (G (1) , G (3) , G (4) ) The figure which shows three examples of the "pattern” corresponding to selection of a subset, and the determination result whether each "pattern" satisfy
  • the figure which shows "the counting up” of the whole graph number (adopted total number) of the graph set represented by ZDD The figure which shows the method of taking out the graph of the arbitrary designated numbers from the graph set represented by ZDD
  • Functional block diagram showing a second embodiment of the present invention Graph showing an example of applying the present invention to a certain gene regulatory network estimation
  • a gene regulatory network applicable to the drug discovery field will be described as an example.
  • the gene regulatory network is expected to be applied to read and understand the mechanism of action of a drug, for example, by expressing the cooperative relationship between genes as a directed graph.
  • FIG. 1 shows a chart showing RNA expression matrix data and a gene regulatory network.
  • RNA expression matrix data is acquired.
  • the RNA expression matrix data may be NGS coverage data or microarray signal data.
  • RNA expression levels of N genes are measured for M cell lines and the like, and data x (m, n) indicates the expression level of gene n in the cell line m. Therefore, the RNA expression matrix D is M ⁇ N numerical matrix data.
  • the relationship between multiple genes can be expressed as a gene regulatory network.
  • the gene regulatory network is represented as a graph G.
  • the graph G is a set of edges (control relationships between nodes (genes) indicated by arrows).
  • g1 ⁇ (A, B), (A, C), (C, D) ⁇
  • g1 ⁇ (A, B), (A, C), (C, D) ⁇
  • the graph G can be estimated from the RNA expression matrix data D for a large number of samples.
  • the graph G should not be a circulant graph. (That is, for example, it should not include a subset such as “(A, B), (B, A)” or “(A, B), (B, C), (C, A)”). If it is expected that the scale-free network property (node order distribution conforms to the power law) will be considered, it may be possible to provide such a restriction.
  • a graph G_1 that seems to be optimal is obtained by some kind of heuristic method for solving the optimization problem, and an evaluation value S_1 is obtained for it.
  • the present invention can be applied to various other combinatorial optimization problems. For example, there is a problem of searching for exclusive and covering gene mutations in cells. Thus, for example, an application for estimating a gene mutation or an action mechanism important for cancer is known.
  • SNP Single Nucleotide Polymorphism
  • Data x (m, n) indicates the presence or absence of SNP mutation n in cell m.
  • SNP refers to a position in DNA (deoxyribonucleic acid) where mutation is likely to occur.
  • a set of loci is represented as a set G.
  • the set G is a set of SNPs.
  • An element of the set G is one of N types of SNPs, and the entire set G corresponds to a combination pattern of SNPs.
  • it is required to exclusively cover M cell lines with the set G. That is, when considering “a set M1 of cells having a mutation in SNP1,” “a set M3 of cells having a mutation in SNP3” and “a set M4 of cells having a mutation in SNP4”, M1, M3, and M4 are mutually overlapping elements. Must be shared, and all cells in M1, M3, and M4 must be covered.
  • greedy hill climbing method for example, greedy hill climbing method, annealing method, tabu search, genetic algorithm, etc. are known, any of which may be used.
  • the present invention can be applied to fields other than so-called bioinformatics.
  • gene regulation network estimation is generalized as a Bayesian network, so it can be used as a technique for measuring various characteristics of a large number of products and converting them into data, and estimating causal relationships between these characteristics.
  • combination optimization problem for example, a knapsack problem and a traveling salesman problem are known and applied in various fields, and the present invention can be applied to any of them.
  • the present invention makes it possible to determine the optimality of the graph G_1 that seems to be optimal.
  • FIG. 2 is a conceptual diagram showing the features of the present invention.
  • an evaluation value evaluation value of the solution in the entire solution space (“score The maximum evaluation value (first maximum evaluation value) Z is also estimated. A specific method for estimating the first maximum evaluation value Z will be described later.
  • the graph G_1 (local solution) considered to be optimal falls within the confidence interval of the estimated first maximum evaluation value Z. If the graph G_1 falls within the confidence interval of the first maximum evaluation value Z, the graph G_1 is the first optimal solution (one of global solutions) in the entire search space (solution space). Can be determined. Such a determination (optimal sufficiency determination as to whether the solution satisfies a sufficient condition) is one of the features of the present invention.
  • the maximum evaluation value (evaluation value) of the evaluation values of the solution on the solution space (subspace) whose distance on the solution space from the graph G_1 is outside a certain range (The second maximum evaluation value) W is estimated.
  • a specific method for estimating the second maximum evaluation value W will be described later.
  • FIG. 3 is a block diagram showing a hardware configuration of the optimum solution determining apparatus according to the present invention.
  • the optimal solution determination apparatus 10 shown in FIG. 3 is configured by a computer, and stores a central processing unit (CPU: Central Processing Unit) 12 that mainly controls the operation of each component and a control program for the device.
  • CPU Central Processing Unit
  • a main memory 14 serving as a work area at the time of execution
  • a graphic board 16 for controlling the display of a monitor device 28 such as a liquid crystal display or a CRT display
  • a communication interface (communication I / F) 18 connected to a network 50
  • Various application software including the optimal solution determination program according to the invention, and the hard disk device 20 for storing the determination result of the optimal solution, which will be described later, the CD-ROM drive 22, and the key operation of the keyboard 30 are detected and input as instructions.
  • Keyboard controller 24 to output to CPU 12 and position
  • the mouse controller 26 detects the state of the mouse 32 as an input device and outputs signals such as the position of the mouse pointer on the monitor device 28 and the state of the mouse 32 to the CPU 12.
  • the network 50 is connected to a database 40 for storing RNA expression matrix data.
  • the RNA expression matrix data is numerical matrix data indicating the RNA expression levels of a plurality of genes (A, B,..., Z) in a plurality of cell lines (samples: X1, X2,..., Xn) as shown in FIG. It is.
  • the RNA expression level is obtained from the sample by NGS (Next Generation Generation Sequencer) (not shown).
  • the optimal solution determination device 10 can access the database 40 via the communication interface 18 and obtain necessary RNA expression matrix data.
  • the RNA expression matrix data is not limited to the one stored in the external database 40, but the RNA expression matrix data is stored in the hard disk device 20, and the RNA expression matrix data stored in the hard disk device 20 is used. You may make it do.
  • FIG. 4 is a functional block diagram showing functions of the CPU 12 of the optimum solution determination apparatus 10 shown in FIG. 3, and is a functional block diagram showing the first embodiment of the present invention.
  • the CPU 12 functions as various processing units by executing the optimal solution determination program stored in the hard disk device 20, and in the first embodiment shown in FIG. 4, the solution extraction unit 100 and the first evaluation value acquisition unit. And the evaluation value acquisition unit 102, the solution input unit 104, the first maximum evaluation value estimation unit 106, the first comparison unit 108, and the first determination unit 110 that function as a second evaluation value acquisition unit. .
  • the solution extraction unit 100 is a part that uniformly extracts a solution (graph) on the solution space of the combinatorial optimization problem (gene regulatory network).
  • a path using ZDD Zero-suppressed binary Decision Diagram
  • the graph G is enumerated and indexed by the enumeration indexing algorithm.
  • ZDD Zero-suppressed binary Decision Diagram
  • the total number of graphs G and arbitrary elements of the set ⁇ G ⁇ can be extracted uniformly.
  • the gene control network consider a graph in which genes are nodes and control relationships are edges.
  • solution here means “acceptable solution (executable solution)”, which is not necessarily an optimal solution, but refers to a solution that is not infeasible. That is, unsuitable solutions (which cannot be executed) are eliminated in advance.
  • the solution extraction unit 100 uses a data structure that uses at least one of “pruning” and “section sharing” to reduce the combinable patterns in the combinatorial optimization problem and enumerate and index them.
  • pruning is to identify whether or not it is determined that a part of the combination is unsuitable without considering the remaining combination among patterns that can be combined in the combination optimization problem, A process for reducing patterns to be identified.
  • section sharing in the frontier method is a pattern to be identified by extracting a common part of a pattern group having a difference in only a part of the combination and combining the remaining combinations. The process to reduce
  • Constraint C is used to determine whether or not “pruning” is “inappropriate”. For example, in the case of a gene control network, if circulation has already occurred due to the employed edge, it is determined that it is inappropriate without considering the remaining edges. Also, in the determination of “whether or not common” for “section sharing”, it is determined whether or not the pattern can be shared in consideration of the considered part of the pattern and the remaining elements under the constraint C. For example, when only the number of edges is considered, if the number of employed edges is the same, “section sharing” can be applied.
  • the “pruning” and “node sharing” algorithms are known as the frontier method in the case of ZDD, and may be followed.
  • random generation may be considered as a means for uniformly extracting solutions. That is, a solution candidate is randomly generated (for example, whether or not there is an edge of the graph is determined by a random number). If the G constraint is not satisfied, regeneration is repeated. If the constraint condition is simple, the constraint may be incorporated at the stage of random generation. For example, when the number of edges is set to a predetermined number or less, an upper limit may be set for the number of edges. On the other hand, for example, when it is desired to prohibit a cyclic graph, it is difficult to apply constraints in simple random generation, so it is possible to determine whether or not a randomly generated graph is circulating.
  • the set partitioning problem is that when a sequence of subsets is given for a certain subset, some of them are selected, and "the selected subsets do not overlap (mutual exclusion)", and "the original subset is exhausted”
  • the problem is whether a pattern (combination) such as “overall coverage” can be created.
  • a set is defined as a “collection of elements”.
  • “G ⁇ " as shown in FIG. 5 shows the whole set, "A, B, C, D" corresponds to the "element”.
  • a “pattern” can be expressed as shown in FIG. In the example of FIG. 7, “patterns” corresponding to selection of three subsets (G (1) , G (3) , G (4) ) are shown.
  • FIG. 8 shows three examples of “patterns” corresponding to selection of a subset, and determination results as to whether or not each “pattern” satisfies a condition.
  • the “pattern” corresponding to the selection of the subset (G (1) , G (3) , G (4) ) satisfies the condition, and the subset (G (2) , G (3) , The “pattern” corresponding to the selection of G (4) ) and subsets (G (3) , G (4) ) does not satisfy the condition.
  • the element “C” is duplicated and does not satisfy the exclusiveness
  • the subset (G (3) , G (4) ) is the element “A”. This is because there is not enough to satisfy the covering property.
  • the combinatorial problem is known to cause a “combined explosion” and it is impossible to search for an optimal solution in a finite time, but due to the reduction techniques of “pruning” and “joint sharing” of the frontier method, An efficient (practical) “enumeration” of a set of combinations that satisfy a predetermined condition can be realized.
  • FIG. 10 is a diagram showing how the combinable patterns are reduced by the “pruning” of the frontier method.
  • FIG. 11 is a diagram showing how the combinable patterns are reduced by “section sharing” of the frontier method.
  • section sharing reduces the number of patterns that can be combined by handling them together if the subsequent patterns are the same.
  • FIG. 12 is a diagram showing a result of reducing combinable patterns by “pruning” and “section sharing” of the frontier method.
  • branches and leaves that should spread to 16 patterns (FIG. 9) can be greatly reduced, and a determination result that completely matches the case of comprehensive determination one by one can be acquired. .
  • FIG. 13 only when passing through “1 branch” indicated by a solid line, it is a graph including the edge (not including “0 branch” indicated by a dotted line or a skipped edge). ”Is adopted, the graph is adopted. (The graph that reaches “0 end” is not adopted.)
  • a graph corresponding to the path of (A ⁇ B ⁇ C) is represented by a route indicated by a thick arrow on the left side of FIG.
  • FIG. 14 is a diagram showing the “counting up” of the total number of graphs (total number of adoption).
  • the total number employed can be calculated by adding “1” to “1 end” indicating the determination result, and counting up from “1 end” to the highest ZDD node in reverse order.
  • the number of grants of each branch point is added from the lower layer ZDD node, and the added number is given to itself. This is repeated up to the highest ZDD node, and the numerical value assigned to the highest ZDD node is the total number of adopted totals (total number of paths reaching “one end”). In the case of this example, the total number of adoption is 20.
  • FIG. 15 is a diagram illustrating a method of extracting a graph having an arbitrary designated number.
  • the graph “No. 12” when the graph “No. 12” is extracted, the graph descends from the highest node according to the thick arrow in FIG. First, the process proceeds from the highest node (A, B) to the branch including the designated number among “0 branch” or “1 branch”.
  • the branch on the “one branch” side is advanced and descends from the highest node (A, B) to the lower node (A, C) (the right node (A, C) in FIG. 14).
  • the branch on the “1 branch” side is advanced, the number on the “0 branch” side is subtracted from the designated number.
  • the number “10” on the “0 branch” side is subtracted from the designated number “12” to become “2”.
  • the solution extraction unit 100 uniformly extracts a solution (graph G) in the solution space of the gene regulatory network by a path enumeration indexing algorithm using ZDD.
  • the total number of graphs G in the solution space can be determined by “counting up”, which is one of the important properties of ZDD, and random numbers (for example, M series ( pseudo-random numbers) using maximal (length) (sequence) are generated, and a graph G corresponding to a designated number designated by each random number is extracted (graph G is uniformly extracted).
  • the evaluation value acquisition unit (first evaluation value acquisition unit) 102 gives an evaluation value to the graph G extracted by the solution extraction unit 100.
  • an evaluation function S (D, G) that quantifies how much the graph G can explain the RNA expression matrix data D is prepared, and the evaluation value acquisition unit 102 corresponds to the extracted graph G.
  • the evaluation value S is acquired based on the evaluation function S (D, G), and the acquired evaluation value S is assigned to the graph G.
  • the evaluation function S (D, G) may be created by the evaluation value acquisition unit 102 based on the RNA expression matrix data stored in the database 40, or may be created based on the RNA expression matrix data in advance, for example, the database 40
  • the evaluation function S (D, G) stored in (1) may be used.
  • the solution acquisition unit 104 acquires a solution (solution candidate) (hereinafter, referred to as “graph G_1”) that seems to be optimal obtained by some heuristic method for solving the optimization problem of the gene regulatory network.
  • the evaluation value S_1 for the input graph G_1 can be acquired by the evaluation value acquisition unit 102.
  • the heuristic method for example, a greedy method, a hill climbing method, an annealing method, a tabu search, a genetic algorithm, and the like are known, and any of them may be used.
  • the heuristic search may be performed by this apparatus or an external apparatus, and the solution acquisition unit 104 acquires a solution candidate (graph G_1) acquired by some heuristic method.
  • the first maximum evaluation value estimation unit 106 is a part that estimates the maximum evaluation value of the solution (graph G) on the solution space.
  • the first maximum evaluation value estimation unit 106 includes a plurality of solutions (first items) uniformly extracted by the solution extraction unit 100.
  • the maximum evaluation value (first maximum evaluation value Z) of the evaluation values when the number exceeding the number of extractions of the first plurality of solutions is assumed is estimated based on the evaluation value of the plurality of one solution). .
  • U and V are natural numbers
  • U ⁇ V solutions (first plurality of graphs G) are uniformly extracted, and an evaluation value S is given to each graph G.
  • U is the block size
  • V is the number of blocks.
  • the first maximum evaluation value estimation unit 106 divides the U ⁇ V graphs G into V blocks, and obtains the segment maximum value among the evaluation values of the U graphs G for each block. Accordingly, it is possible to obtain V maximum segment values. Then, the maximum evaluation value (first maximum evaluation value Z) is estimated with maximum likelihood, assuming that the V segment maximum values follow a generalized extreme value distribution (GEV).
  • GEV generalized extreme value distribution
  • the first maximum evaluation value is accompanied by statistical support.
  • the graph ⁇ G ⁇ on the solution space is originally a finite set and degenerates strictly. However, since the total number of graphs G is sufficiently large, continuous distribution approximation can be applied. In that case, since there is a clear upper limit in the evaluation value S, it is expected that the evaluation value S becomes a Gumbel type by setting appropriate U and V, and the true first maximum evaluation value Z can be estimated with a confidence interval.
  • the first comparison unit 108 acquires the evaluation value S_1 corresponding to the solution candidate (graph G_1) acquired by the solution acquisition unit 104 from the evaluation value acquisition unit 102, and has the confidence interval with the estimated evaluation value S_1.
  • the maximum evaluation value Z of 1 is compared.
  • the first determination unit 110 determines whether or not the evaluation value S_1 of the graph G_1 falls within the confidence interval of the first maximum evaluation value Z based on the comparison result by the first comparison unit 108. If it falls within the confidence interval of the first maximum evaluation value Z, it can be seen that the graph G_1 is one of the first optimal solutions (“graph G_1 is sufficient”) in the entire solution space.
  • FIG. 16 is a functional block diagram showing functions of the CPU 12 of the optimum solution determining apparatus 10 shown in FIG. 3, and is a functional block diagram showing a second embodiment of the present invention.
  • the same reference numerals are given to the portions common to the first embodiment shown in FIG. 4, and detailed description thereof is omitted.
  • a second maximum evaluation value estimation unit 112 a second comparison unit 114, and a second determination unit 116 are mainly added as compared with the first embodiment.
  • the second embodiment makes it possible to determine that there is no other solution candidate that is equally likely other than the solution candidate (G_1).
  • the solution extraction unit 100 selects a solution whose distance in the solution space from the graph G_1 is outside a certain range. Construct ZDD for enumeration indexing. This ZDD construction can be realized by counting the common or non-common edges of the graph G and the graph G_1 during construction in the frontier method for constructing the ZDD.
  • the solution extraction unit 100 is a second plurality of solutions (graph G) in the solution space of the gene control network by the path enumeration indexing algorithm using the reconstructed ZDD, and the distance in the solution space from the graph G_1 A second plurality of graphs G outside a certain range are uniformly extracted.
  • the second maximum evaluation value estimation unit 112 evaluates when the number exceeding the extracted number of the second plurality of solutions is assumed based on the evaluation values respectively corresponding to the second plurality of uniformly extracted graphs.
  • the maximum evaluation value (second maximum evaluation value W) of the values is estimated.
  • the second maximum evaluation value W is estimated by a method similar to the estimation of the first maximum evaluation value Z by the first maximum evaluation value estimation unit 106. That is, if P and Q are natural numbers, P ⁇ Q solutions (second plurality of graphs G) are uniformly extracted, and an evaluation value S is given to each graph G.
  • P is the block size and Q is the number of blocks.
  • P and Q may be the same as U and V uniformly extracted when the first maximum evaluation value Z is estimated, or may be different.
  • the second maximum evaluation value estimation unit 112 divides the P ⁇ Q graphs G into Q blocks, and acquires the segment maximum value among the evaluation values of the U graphs G for each block. Therefore, it is possible to obtain Q pieces of maximum classification values. Then, the second maximum evaluation value W is estimated with maximum likelihood, assuming that the Q segment maximum values follow the general extreme value distribution.
  • the second comparison unit 114 acquires the evaluation value S_1 corresponding to the solution candidate (graph G_1) acquired by the solution acquisition unit 104 from the evaluation value acquisition unit (second evaluation value acquisition unit) 102, and the acquired evaluation value S_1 is compared with the estimated second maximum evaluation value W with the confidence interval.
  • the second determination unit 116 determines that the evaluation value S_1 of the graph G_1 based on the comparison result by the second comparison unit 114 is the second maximum evaluation value W (the range of the second maximum evaluation value W with a confidence interval). It is determined whether or not it exceeds.
  • the first maximum evaluation value Z is an estimate of the maximum value of the entire solution space
  • the second maximum evaluation value W is an estimate of the maximum value of the subspace
  • W ⁇ Z is Obviously (depending on the sample size or the like, there may be a case where W> Z stochastically).
  • the graph G_1 is a graph that gives the first maximum evaluation value, and in the range away from the graph G_1 in the solution space, it is equal to or greater than It can be determined that there is no graph structure that can explain the RNA expression matrix data D.
  • the distance within a certain range to be separated from the graph G_1 needs to be set in advance, but this may be set from the characteristics of the graph, may be set empirically, for example, S_1 >> W You may increase the distance which gradually separates to a certain point. For example, an appropriate distance may be repeatedly searched efficiently from a sufficiently large distance setting value by a two-minute search or the like.
  • the graph G_1 is the only optimal solution (“graph G_1 is necessary and sufficient”) in the entire solution space. That is, it can be determined that there is no graph corresponding to the optimum value other than the graph G_1 obtained by the heuristic search.
  • the solution candidate in a range away from the graph G may be used as a means for directly comparing the graph G with the solution candidate.
  • the present invention is a method supported by a statistical basis.
  • FIG. 17 is a graph showing an example in which the present invention is applied to a certain gene regulatory network estimation.
  • the horizontal axis of the graph shown in FIG. 17 is the degree of divergence (distance), and the vertical axis is the evaluation value.
  • a dotted line is a reaching evaluation value (first maximum evaluation value Z) for the acquired graph.
  • the white line is an estimated value of the optimum value for each degree of divergence.
  • the “degree of deviation” here is generally a distance or an index according to it.
  • the Hamming distance may be used.
  • other distance indices such as edit distance, or indices such as quasi-distance and half-distance may be used.
  • the achievement evaluation value is determined to be the optimum value.
  • the value corresponding to the estimated range gradually decreases, and the estimated range (second maximum evaluation value W) deviates from the reached evaluation value at the divergence degree 5 to 6, so the divergence degree 5 to 6 It was determined that there was no graph having an evaluation value equal to or higher than the acquired graph in the above range.
  • the solid line indicates the actual measurement value corresponding to the solution candidate obtained by the heuristic search, and indicates that the present invention can correctly estimate the actual situation.
  • FIG. 18 is a flowchart showing the first embodiment of the optimal solution determination method according to the present invention.
  • the solution extraction unit 100 shown in FIG. 4 uniformly extracts a solution (graph G) on the solution space of the gene regulatory network that is one of the combinatorial optimization problems (step S10 (first step)).
  • the graph G is uniformly extracted by the path enumeration indexing algorithm using ZDD as described above.
  • the evaluation value acquisition unit 102 assigns an evaluation value S to the uniformly extracted graph G (step S12).
  • an evaluation function S (D, G) that quantifies how much the graph G can explain the RNA expression matrix data D is prepared, and the evaluation value acquisition unit 102 applies the uniformly extracted graph G to the graph G.
  • the corresponding evaluation value S is acquired based on the evaluation function S (D, G), and the acquired evaluation value S is assigned to the graph G.
  • the first maximum evaluation value estimation unit 106 acquires the first plurality of evaluation values S given in step S12 for the plurality of uniformly extracted graphs G (first plurality of solutions) (Ste S14, second step), based on the acquired first plurality of evaluation values S, the maximum evaluation value (evaluation value among the evaluation values when the number exceeding the extracted number of the first plurality of graphs G is assumed ( The first maximum evaluation value) Z is estimated (step S16, third step).
  • the first plurality of graphs G are U ⁇ V graphs, where U and V are natural numbers, respectively, and the first maximum evaluation value estimation unit 106 includes the U ⁇ V graphs G. Is divided into V blocks, V pieces of the maximum value of the evaluation values of the U graphs G for each block are obtained, and the maximum number of pieces follows the general extreme value distribution using the V pieces of maximum values.
  • the maximum evaluation value (first maximum evaluation value Z with confidence classification) is estimated with maximum likelihood.
  • the solution acquisition unit 104 acquires an optimal solution candidate (graph G_1) obtained by some heuristic method for solving the optimization problem of the gene regulatory network (step S18, fourth step).
  • the heuristic method for example, a greedy method, a hill climbing method, an annealing method, a tabu search, a genetic algorithm, and the like are known, and any of them may be used.
  • Evaluation value acquisition unit 102 acquires evaluation value S_1 for the acquired graph G_1 (step S20, fifth step).
  • the first comparison unit 108 compares the evaluation value S_1 acquired in step S20 with the first maximum evaluation value Z with a confidence interval estimated in step S16 (step S22, sixth step).
  • the first determination unit 110 determines whether or not the evaluation value S_1 of the graph G_1 falls within the range of the first maximum evaluation value Z with a confidence interval based on the comparison result by the first comparison unit 108 ( Step S24, seventh step). That is, when the evaluation value S_1 falls within the confidence interval of the first maximum evaluation value with a confidence interval, the first determination unit 110 satisfies the sufficient condition as an optimal solution ( It is determined that there is an optimality.
  • the determination result by the first determination unit 110 is displayed on the monitor device 28 shown in FIG. 3, stored in the hard disk device 20, or printed out to a printer (not shown) (step S26, eighth step).
  • the determination result by the first determination unit 110 is not limited to the presence or absence of optimal sufficiency. If there is no optimal sufficiency, the difference between the evaluation value S_1 and the first maximum evaluation value Z is calculated on the solution space. It may be converted into a distance, and it may be determined how far the currently estimated graph G_1 is from the true optimum solution (the solution corresponding to the first maximum evaluation value Z).
  • the solution candidate (graph G_1) searched by the heuristic search has optimal sufficiency as a true optimal value.
  • the solution candidate may be used in consideration of the fact that there is a separate optimum solution. That is, the information on the optimal sufficiency determination is useful regardless of the success or failure of the optimal sufficiency determination.
  • a heuristic search method there are a first search method with low calculation cost (short search time) but low solution accuracy, and a high calculation cost (long search time) than the first search method.
  • the second search method with high accuracy of the first search method is prepared, and the optimal sufficiency of the first solution candidate (graph G_1) first searched by the first search method is determined, and the determination of the optimal sufficiency fails Only when it is determined, the optimal sufficiency of the second solution candidate (graph G_1) searched by the second search method is determined.
  • search methods may be switched between heuristic search methods, or may be switched to an approximation algorithm with a certain degree of guarantee of approximation or a method for obtaining an exact solution.
  • Three or more search methods may be prepared and sequentially switched. Switching of search methods is not limited to the method itself, and may be realized by convergence determination of the same search method. For example, in a search method that increases the accuracy by repetitive search, the search result of a predetermined number of times may be determined by the first maximum evaluation value Z, and the search may be repeated until it reaches the confidence interval of the first maximum evaluation value Z.
  • the first maximum evaluation value Z for determining the optimum sufficiency may be acquired in advance.
  • FIG. 19 is a flowchart showing the second embodiment of the optimum solution judging method according to the present invention, and in particular, the optimum necessity judgment performed after the optimum sufficiency judgment according to the first embodiment shown in FIG. It shows about processing.
  • the solution extraction unit 100 (second solution extraction unit) shown in FIG. 16 is a second plurality of graphs G on the solution space of the gene regulatory network, and the graph that has succeeded in determining the optimal sufficiency
  • a plurality of (second plurality) graphs G whose distance in the solution space from G_1 is outside a certain range are uniformly extracted (step S30, ninth step).
  • the uniform extraction of the second plurality of graphs G is performed by constructing a ZDD that enumerates and indexes solutions whose distance in the solution space from the graph G_1 is outside a certain range, and performs path enumeration indexing using the constructed ZDD. This can be done with an algorithm.
  • the evaluation value acquisition unit 102 acquires evaluation values S for the second plurality of uniformly extracted graphs G (step S32, tenth step). Acquisition of the evaluation value S for the graph G can be performed in the same manner as step S12 of the first embodiment shown in FIG.
  • the second maximum evaluation value estimation unit 112 evaluates when the number exceeding the number of extracted second graphs G is assumed based on the second plurality of evaluation values S acquired in step S32.
  • the maximum evaluation value (second maximum evaluation value) W of the values is estimated (step S34, eleventh step).
  • the second plurality of graphs G are P ⁇ Q graphs, where P and Q are natural numbers, respectively, and the second maximum evaluation value estimation unit 112 has P ⁇ Q graphs G. Is divided into Q blocks, Q is obtained as the maximum value of the evaluation values of the P graphs G for each block, and the maximum value of the division follows the general extreme value distribution using the Q maximum values of the division.
  • the maximum evaluation value (second maximum evaluation value W with confidence classification) is estimated to the maximum likelihood.
  • the second comparison unit 114 compares the evaluation value S_1 of the graph G_1 that has succeeded in determining the optimal sufficiency with the second maximum evaluation value W with the confidence interval estimated in Step S34 (Step S36, 12th step).
  • the second determination unit 116 determines whether or not the evaluation value S_1 of the graph G_1 exceeds the range of the second maximum evaluation value W with a confidence interval based on the comparison result by the second comparison unit 114. (Step S38, 13th step). That is, when the evaluation value S_1 exceeds the range of the second maximum evaluation value with the confidence interval, the second determination unit 116 has no other solution equivalent to the graph G_1, and the graph G_1 It is determined that the necessary condition as the only optimal solution is satisfied (there is an optimal necessity).
  • the determination result by the second determination unit 116 is displayed on the monitor device 28 shown in FIG. 3, stored in the hard disk device 20, or printed out to a printer (not shown) (step S40, 14th step).
  • the search is performed to determine the optimal necessity of the graph G_1. It can be confirmed that there is no solution having an equivalent evaluation value in the solution space away from the graph G_1.
  • the third embodiment of the optimal solution determination method according to the present invention includes processing when the optimal sufficiency determination fails in the second embodiment shown in FIG.
  • step S30 (fifteenth step) shown in FIG. 19, a certain range from the graph G_1 is enlarged on the solution space, and on the solution space outside the enlarged certain range.
  • a third plurality of graphs G_3 (third plurality of solutions) are uniformly extracted.
  • step S32, 16th step acquisition of the third plurality of evaluation values corresponding to the third plurality of graphs G extracted in the solution space outside the fixed range
  • step S34, 17th step estimation of the third maximum evaluation value W when assuming the number of solutions exceeding the number of the third plurality of graphs G
  • step S_1 of the graph G_1 and the third Are compared with the maximum evaluation value W (step S36), the determination of the optimal necessity of the graph G_1 (step S38, the 18th step), etc.
  • the above-described processing for expanding a certain range may be repeated multiple times while gradually expanding the certain range until the optimal sufficiency determination is successful.
  • the fourth embodiment of the optimal solution determination method according to the present invention includes other processing when the optimal sufficiency determination fails in the second embodiment shown in FIG.
  • the graph G_2 is obtained by performing a heuristic search for a solution candidate (graph G_2) in a solution space whose distance from the graph G_1 is outside a certain range, and the evaluation value S_2 of the graph G_2 is equivalent to W ⁇ Z This is done by obtaining G_2.
  • step S30 a fourth plurality of solutions whose distances in the solution space from the graph G_1 and the graph G_2 are outside a certain range are uniformly extracted (steps S44, 20th). Step).
  • step S32 shown in FIG. 19 to obtain a fourth plurality of evaluation values corresponding to the fourth plurality of graphs G extracted on the solution space outside the fixed range expanded as described above ( (Step S32, 21st step), estimation of the fourth maximum evaluation value W when assuming the number of solutions exceeding the number of the fourth plurality of graphs G based on the fourth plurality of evaluation values (Step S34) 22nd step), comparing the evaluation value S_1 of the graph G_1 with the second maximum evaluation value W (step S36), determining the optimum necessity of the graph G_1 (step S38, 23rd step), etc. Execute.
  • a certain range is determined by the request under the combinatorial optimization problem.
  • the solution cannot be extracted uniformly from the entire solution space. This is a case where uniform extraction means cannot be secured in the first place, such as a case where ZDD of the entire solution space cannot be constructed.
  • FIG. 21 is a flowchart showing the fifth embodiment of the optimum solution determination method according to the present invention.
  • the solution space of the combinatorial optimization problem includes the first solution space in the above-described constraint condition C (first constraint condition) and the second solution space in the constraint condition C + (second constraint condition). .
  • the first solution candidate (graph G_1) is acquired from the solutions belonging to the first solution space ( ⁇ G ⁇ in the constraint condition C) (step S50, step 24).
  • the graph G_1 can be searched by a heuristic method or the like.
  • the solutions in the second solution space (graph ⁇ G + ⁇ in the constraint condition C + ) are listed, and the neighborhood solution G_1 + in the constraint condition C + of the graph G_1 acquired in step S50 is searched (step S52, 28th step).
  • the ZDD in the constraint condition C + can be described, the ZDD is constructed by adding “the distance in the solution space from G_1 is within a certain range” to the constraint condition C + and constructing the ZDD of the minimum distance.
  • the neighborhood solution G_1 + can be searched by selecting one.
  • evaluation value S_1 + of the neighborhood solution G_1 + is acquired (step S54, 29th step).
  • Evaluation value S_1 + acquisition for neighborhood solutions G_1 + can be carried out in the same manner as step S12 of the first embodiment shown in FIG. 18.
  • Step S56 a neighborhood solutions G_1 + evaluation value S_1 + acquired in S54, is compared with the maximum evaluation value Z of the fifth with confidence intervals (step S56).
  • the fifth maximum evaluation value Z can be estimated by the above-described method based on the uniform extraction of the graph ⁇ G + ⁇ in the constraint condition C + in the second solution space. That is, a plurality of solutions on the second solution space are uniformly extracted as a fifth plurality of solutions (graph ⁇ G + ⁇ in the constraint C + ) (25th step), and the fifth plurality of solutions Fifth plurality of evaluation values corresponding to each are acquired (26th step), and the number of solutions exceeding the number of fifth plurality of solutions is assumed based on the acquired fifth plurality of evaluation values. In this case, the fifth maximum evaluation value is estimated (27th step).
  • the optimal sufficiency of the graph G_1 is determined (step S58, the 30th step). . That is, it is determined whether or not the evaluation value S_1 + falls within the confidence interval of the first maximum evaluation value Z. The evaluation value S_1 + falls within the confidence interval of the first maximum evaluation value Z, and the neighborhood solution G_1 + If it is determined that there is optimal sufficient, the graph G_1 is also determined to have optimal sufficient. This is because, if the neighborhood solution G_1 + is optimal, it can be ensured that at least “the closest solution is found from the optimal solution in the case where more restrictive conditions are imposed”.
  • FIG. 22 is a diagram schematically illustrating the first solution space under the constraint condition C and the second solution space under the constraint condition C + .
  • the second solution space is “C + ⁇ C” because the constraint is added (FIG. 22). This is because the constraint included in the constraint condition C is a precondition of the model or a reasonable constraint based on prior knowledge, and thus should not be basically removed. Needless to say, it is desirable that the added constraint is assumed to have a certain degree of validity with respect to the problem by prior knowledge or the like (although it cannot be guaranteed even if the complete validity is not guaranteed).
  • the difference between the constraints C and C + and the minimum distance d between the graph G_1 and the neighborhood solution G_1 + indicate the degree of divergence between the substitute neighborhood solution G_1 + and the actually used solution.
  • the search for the neighborhood solution G_1 + if it can be said that the neighborhood of the minimum distance d of the reaching solution is searched by, for example, a heuristic method, it is considered that there is no practical problem.
  • the influence of the difference between the constraint conditions C and C + depends on how much the evaluation value changes accordingly.
  • the reliability determination for the optimality determination of the graph G_1 is high by substituting the optimality determination of the neighborhood solution G_1 + You may think. From the above, the effect of the present invention can be exerted with respect to a wider combination optimization problem.
  • the optimal solution determination apparatus 10 of this embodiment is merely an example, and the present invention can be applied to other configurations.
  • Each functional configuration can be appropriately realized by arbitrary hardware, software, or a combination of both.
  • an optimal solution determination program that causes a computer to execute the processing in each unit of the above-described optimal solution determination device 10 and a computer-readable recording medium (non-temporary recording medium) that records such an optimal solution determination program
  • the present invention can be applied.
  • the solution extraction unit 100 the evaluation value acquisition unit 102, the solution acquisition unit 104, the first maximum evaluation value estimation unit 106, the first comparison unit 108, the first determination unit 110, and the like.
  • the hardware structure of the processing unit (processing unit) that executes the various types of processing is the following types of processors.
  • the circuit configuration can be changed after manufacturing a CPU (Central Processing Unit) or FPGA (Field Programmable Gate Array) that is a general-purpose processor that functions as various processing units by executing software (programs).
  • a CPU Central Processing Unit
  • FPGA Field Programmable Gate Array
  • dedicated logic circuits such as programmable logic devices (Programmable Logic Devices: PLDs) and ASICs (Application Specific Specific Integrated Circuits) that have specially designed circuit configurations to execute specific processing. It is.
  • One processing unit may be configured by one of these various processors, or may be configured by two or more processors of the same type or different types (for example, a plurality of FPGAs or a combination of CPU and FPGA). May be. Further, the plurality of processing units may be configured by one processor. As an example of configuring a plurality of processing units with one processor, first, as represented by a computer such as a client or server, one processor is configured with a combination of one or more CPUs and software. There is a form in which the processor functions as a plurality of processing units.
  • SoC system-on-chip
  • a form of using a processor that realizes the functions of the entire system including a plurality of processing units with a single IC (integrated circuit) chip. is there.
  • various processing units are configured using one or more of the various processors as a hardware structure.
  • circuitry circuitry in which circuit elements such as semiconductor elements are combined.
  • the present invention is also an optimum solution determination apparatus having a processor, wherein the processor uniformly extracts a plurality of solutions on the solution space of the combinatorial optimization problem as a first plurality of solutions and uniformly extracts the first solution.
  • the first plurality of evaluation values corresponding to each of the plurality of one solution are obtained, and the number of solutions exceeding the number of the first plurality of solutions is assumed based on the obtained first plurality of evaluation values.
  • the evaluation value corresponding to the acquired solution candidate is the first maximum evaluation value.
  • the optimal solution determination apparatus which determines whether it falls in within the confidence interval of is included.

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Health & Medical Sciences (AREA)
  • General Physics & Mathematics (AREA)
  • Mathematical Physics (AREA)
  • Data Mining & Analysis (AREA)
  • Biophysics (AREA)
  • Bioinformatics & Cheminformatics (AREA)
  • Bioinformatics & Computational Biology (AREA)
  • Evolutionary Biology (AREA)
  • Software Systems (AREA)
  • General Engineering & Computer Science (AREA)
  • General Health & Medical Sciences (AREA)
  • Medical Informatics (AREA)
  • Evolutionary Computation (AREA)
  • Artificial Intelligence (AREA)
  • Mathematical Analysis (AREA)
  • Computational Mathematics (AREA)
  • Biotechnology (AREA)
  • Pure & Applied Mathematics (AREA)
  • Mathematical Optimization (AREA)
  • Spectroscopy & Molecular Physics (AREA)
  • Computational Linguistics (AREA)
  • Computing Systems (AREA)
  • Genetics & Genomics (AREA)
  • Molecular Biology (AREA)
  • Operations Research (AREA)
  • Databases & Information Systems (AREA)
  • Business, Economics & Management (AREA)
  • Physiology (AREA)
  • Algebra (AREA)
  • Strategic Management (AREA)
  • Human Resources & Organizations (AREA)
  • Economics (AREA)
  • Proteomics, Peptides & Aminoacids (AREA)
  • Chemical & Material Sciences (AREA)
  • Biomedical Technology (AREA)
  • Analytical Chemistry (AREA)

Abstract

組合せ最適化問題における解の最適性をコンピュータにより判定する最適解判定方法であって、組合せ最適化問題の解空間上の複数の解を第１の複数の解として一様抽出し、一様抽出した第１の複数の解のそれぞれに対応する第１の複数の評価値に基づいて、第１の複数の解の個数を超える個数の解を想定した場合の最大評価値を第１の最大評価値Ｚとして推定する。そして、解空間に属する解候補（グラフＧ_1）を入力すると（ステップＳ１８）、そのグラフＧ_1に対応する評価値Ｓ_1を取得し、取得した評価値Ｓ_1と第１の最大評価値Ｚとを比較し、入力したグラフＧ_1の評価値Ｓ_1が第１の最大評価値Ｚの信頼区間内に入るか否か（グラフＧ_1が第１の最適解か否か）を判定する。

Description

最適解判定方法、最適解判定プログラム及び最適解判定装置
　本発明は最適解判定方法、最適解判定プログラム及び最適解判定装置に係り、特に組合せ最適化問題における解の最適性を判定する技術に関する。
　近年、ビッグデータに対するデータマイニングなど、データ解析のニーズが高まってきている。重要なデータ解析分野のひとつとして、組合せ最適化問題がある（例えば古くから知られている巡回セールスマン問題などが含まれる）。
　組合せ最適化問題は、ＮＰ（Non-deterministic Polynomial time）完全あるいはＮＰ困難な難しい問題が多く含まれる。即ち、一般的に問題の規模が大きくなると計算量が指数関数以上のオーダーで爆発するので、網羅的な全探索による解決がほとんど不可能である。
　そこで、ヒューリスティックに近似解を求める手法が多く開発され、適用されている。しかし、全探索をしていないので、得られた解が「局所解」に留まり、真の最適解かどうかを判定できないという問題がある。
　この問題に直接応えるため、解に何らかのスコアを付与し、統計解析を併用することで、解の最適性を見極めようという研究があり、例えば、特に極値統計解析によって最大スコアを推定する手法などがある（非特許文献１）。
　しかしながら、場合によっては複雑な制約条件が課せられる組合せ最適化問題においては、良質な統計解析の前提となる解空間から解候補を大量に一様抽出する手段の提供が難しく、実用化には至っていないと考えられる（非特許文献２）。
　また、組合せ最適化問題に関連して、ＺＤＤ（Zero-suppressed Binary Decision Diagram）と呼ばれるデータ構造があり、フロンティア法と称される構築アルゴリズムによって、非常に大規模な組合せ集合を効率的に列挙索引化できることが分かり、近年盛んに研究されている（非特許文献３）。ＺＤＤでは、組合せ最適化問題の解候補全体を描出し、そこから効率的に解候補を一様抽出できる。
　また、創薬分野では、近年発達したＮＧＳ（Next Generation Sequencer：次世代シーケンサ）などにより、大量の遺伝子データ、例えばＲＮＡ（ribonucleic acid）発現行列データが取得できるようになった。そうして得られたビッグデータの解析は、バイオインフォマティクスとして注目されている。例えば、生体機能を踏まえた薬剤の作用機序などを解明しようという試みがある。そのひとつに、遺伝子制御ネットワークの推定がある。遺伝子制御ネットワークとは、遺伝子が相互に発現量を調節するシステムを、ベイジアンネットワークなどの確率的グラフモデルとして捉える解析手法である。
　特許文献１には、遺伝子発現の時系列研究から遺伝子間のネットワーク関係を推定するためのノンパラメトリック回帰によるベイジアンモデルを使用する技術が開示されている。
国際公開第２００４／０４７０２０号
Golden, B.L. and Alt, F.B. "Interval estimation of a global optimum for large combinatorial problems", Naval Research Logistics Quarterly, 26, 69-77 (1979) Giddings, A.P., Rardin, R.L., and Uzsoy, R. "Statistical optimum estimation techniques for combinatorial optimization problems: a review and critique", Journal of Heuristics, 20(3), 329-358 (2014) Iwashita, H., Kawahara, J., and Minato, S. "ZDD-Based Computation of the Number of Paths in a Graph", TCS Technical Report, TCS-TR-A-12-60, Hokkaido University, September 18, 2012.
　遺伝子制御ネットワークでは、遺伝子をノード、制御関係をエッジとするグラフを考える。そして、グラフ構造によって得られているＲＮＡ発現行列データをどの程度説明し得るかを計算して、得られているデータに最も適合するグラフを探索しようという、グラフ探索（グラフマイニング）問題を解きたい。しかしながら、遺伝子数Ｎに対して、可能なグラフ数は2^(N^2)通りあり、Ｎの増大につれて超指数関数的に発散する。また、ベイジアンネットワークモデルでは、非循環有向グラフ制約（DAG制約：Directed acyclic graph制約）という複雑な制約条件も考えなければならない。遺伝子数がある程度大きくなると、ＺＤＤでも解空間全体を描出するには困難を伴う。
　しかも、創薬プロジェクトは、人体に供される可能性がある薬剤を生み出そうというもので、結果の判明まで長い時間と莫大なコストを要し、産業上特に、解析結果の妥当性に対する関心が高い。よって、取得したＲＮＡ発現行列データを巧く説明するグラフ構造の最適性が判定できれば、産業上非常に有用である。
　特許文献１に記載の遺伝子制御ネットワークを構築する方法は、いくつかの遺伝子に対して遺伝子発現の経時変化データを取得し、ベイジアン推定法の修正と、発現遺伝子間の原因及び結果の関係を判断する。ベイジアン推定法の修正は、経時変化データを用いて発現遺伝子間の因果関係を推定することを含む。特許文献１に記載の方法は、発現の遅い遺伝子の変化は、発現の早い遺伝子の変化の原因となる可能性があまりないという仮定に基づいてベイジアン推定及びノンパラメトリック回帰を修正して信頼性の高いネットワーク解を提供するものであるが、上記仮定が全ての遺伝子制御ネットワークに適用できる保証はない。
　本発明はこのような事情に鑑みてなされたもので、組合せ最適化問題における解の最適性の判定を効率的かつ精度よく行うことができる最適解判定方法、最適解判定プログラム及び最適解判定装置を提供することを目的とする。
　上記目的を達成するために一の態様に係る発明は、組合せ最適化問題における解の最適性をコンピュータにより判定する最適解判定方法であって、組合せ最適化問題における解の最適性をコンピュータにより判定する最適解判定方法であって、組合せ最適化問題の解空間上の複数の解を第１の複数の解として一様抽出する第１のステップと、第１のステップにより一様抽出した第１の複数の解のそれぞれに対応する第１の複数の評価値を取得する第２のステップと、取得した第１の複数の評価値に基づいて、第１の複数の解の個数を超える個数の解を想定した場合の最大評価値を推定し、第１の最大評価値とする第３のステップと、解空間に属する解のうち少なくとも１つの解を解候補として取得する第４のステップと、解候補に対応する評価値を取得する第５のステップと、第５のステップにより取得した解候補に対応する評価値が第１の最大評価値の信頼区間内に入るか否かを判定する第６のステップと、第６のステップにおいて解候補に対応する評価値が第１の最大評価値の信頼区間内に入ると判定された場合に、解候補を第１の最適解とする第７のステップと、を含む。
　本発明の一の態様によれば、組合せ最適化問題の解空間上の複数の解を第１の複数の解として一様抽出し、一様抽出した第１の複数の解の個数を超える個数の解を想定した場合の最大評価値を第１の最大評価値として推定する。そして、解空間に属する解のうち少なくとも１つの解を解候補として取得すると、その解候補に対応する評価値を取得し、取得した評価値が第１の最大評価値の信頼区間内に入るか否か（解候補が第１の最適解の１つか否か）を判定する。これにより、組合せ最適化問題の解の統計的な最適性判定が可能になり、かつ解が十分条件を満たすか否かの最適性判定が可能である。
　本発明の他の態様に係る最適解判定方法において、第１の複数の解は、Ｕ，Ｖをそれぞれ自然数とすると、Ｕ×Ｖ個の解であり、第３のステップは、Ｕ×Ｖ個の解をＶ個のブロックに分け、ブロック毎にＵ個の解の評価値の区分最大値をＶ個取得し、Ｖ個の区分最大値を用いて、区分最大値が一般極値分布に従うものとして第１の最大評価値を推定することが好ましい。ブロック毎のＵ個の解の評価値の区分最大値をＶ個取得し、Ｖ個の区分最大値を用いて、区分最大値が一般極値分布（ＧＥＶ：generalized extreme value distribution）に従うものとして最大評価値（第１の最大評価値）を最尤推定する。
　本発明の更に他の態様に係る最適解判定方法において、第７のステップによる判定結果を出力する第８のステップを更に含むことが好ましい。
　本発明の更に他の態様に係る最適性判定方法において、演算コストは小さいが解の精度が低い第１の探索法と、第１の探索法よりも演算コストは大きいが解の精度が高い第２の探索法とを有し、第４のステップは、最初に第１の探索法により探索された第１の解候補を入力し、第１の解候補の評価値が第１の最大評価値の信頼区間内に入らない場合のみ、第２の探索法により探索された第２の解候補を入力することが好ましい。
　最初に演算コストは小さいが解の精度が低い第１の探索法により探索された第１の解候補を入力し、その第１の解候補が十分条件を満たさず十分性判定に失敗した場合、第１の探索法によるヒューリスティック探索が不十分なために、他の最適解が存在することが示唆されたことになる。その場合、演算コストは大きいが解の精度が高い第２の探索法による第２の解候補の探索に切り替え、第２の探索法により探索される第２の解候補の十分性を判定する。
　本発明の更に他の態様に係る最適解判定方法において、解空間上の解であって、第１の最適解からの解空間上の距離が一定範囲外の複数の解を第２の複数の解として一様抽出する第９のステップと、第２の複数の解のそれぞれに対応する第２の複数の評価値を取得する第１０のステップと、第２の複数の評価値に基づいて、第２の複数の解の個数を超える個数の解を想定した場合の最大評価値を推定し、第２の最大評価値とする第１１のステップと、第１の最適解の評価値を取得する第１２のステップと、第１の最適解の評価値が、第２の最大評価値を超えているか否かを判定する第１３のステップと、を含むことが好ましい。
　本発明の更に他の態様によれば、第１の最適解からの解空間上の距離が一定範囲外の複数の解を第２の複数の解として一様抽出し、一様抽出した第２の複数の解の個数を超える個数の解を想定した場合の最大評価値を第２の最大評価値として推定する。そして、第１の最適解の評価値が第２の最大評価値を超えているか否かを判定する。これにより、組合せ最適化問題の解が必要条件を満たすか否かの最適性判定が可能であり、第１の最適解の評価値が第２の最大評価値を超えている場合には、第１の最適解は必要十分条件を満たし、解空間上で最適な解であり、同等解も他には存在しないことになる。
　本発明の更に他の態様に係る最適解判定方法において、第２の複数の解は、Ｐ，Ｑをそれぞれ自然数とすると、Ｐ×Ｑ個の解であり、第１１のステップは、Ｐ×Ｑ個の解をＱ個のブロックに分け、ブロック毎のＰ個の解の評価値の区分最大値をＱ個取得し、Ｑ個の区分最大値を用いて、区分最大値が一般極値分布に従うものとして第２の最大評価値を推定することが好ましい。ブロック毎のＰ個の解の評価値の区分最大値をＱ個取得し、Ｑ個の区分最大値を用いて、区分最大値が一般極値分布に従うものとして最大評価値（第２の最大評価値）を最尤推定する。
　本発明の更に他の態様に係る最適解判定方法において、第１３のステップによる判定結果を出力する第１４のステップを更に含むことが好ましい。
　本発明の更に他の態様に係る最適解判定方法において、第１３のステップにより第１の最適解の評価値が第２の最大評価値を超えていないと判定されると、一定範囲を拡大し、拡大した一定範囲外の複数の解を第３の複数の解として一様抽出する第１５のステップと、第３の複数の解のそれぞれに対応する第３の複数の評価値を取得する第１６のステップと、第３の複数の評価値に基づいて、第３の複数の解の個数を超える個数の解を想定した場合の最大評価値を推定し、第３の最大評価値とする第１７のステップと、第１の最適解の評価値が、第３の最大評価値を超えているか否かを判定する第１８のステップと、を更に含むことが好ましい。
　最適解の１つとして探索された第１の最適解の評価値が第２の最大評価値を超えていない場合、同等の最適解が存在する可能性があるが、この場合には、第１の最適解から離れた解空間（拡大した一定範囲外の解空間）から一様抽出した第３の複数の解に基づいて第３の最大評価値を再度推定する。これにより、探索された最適解から離れた解空間には、同等の最適解が存在しないことを確認することができる。
　本発明の更に他の態様に係る最適解判定方法において、第１３のステップにより第１の最適解の評価値が第２の最大評価値を超えていないと判定されると、解空間上の解であって、かつ第１の最適解から一定距離離れた第４の解候補を取得する第１９のステップと、解空間上において、第１の最適解及び第４の解候補のそれぞれからの解空間上の距離が一定範囲外の第４の複数の解を一様抽出する第２０のステップと、第４の複数の解のそれぞれに対応する第４の複数の評価値を取得する第２１のステップと、第４の複数の評価値に基づいて、第４の複数の解の個数を超える個数の解を想定した場合の最大評価値を推定し、第４の最大評価値とする第２２のステップと、第１の最適解の評価値が、第４の最大評価値を超えているか否かを判定する第２３のステップと、を更に含むことが好ましい。
　最適解の１つとして探索された第１の最適解の評価値が第２の最大評価値を超えていない場合、第１の最適解と、第１の最適解から一定距離離れた第４の解候補との両方から、それぞれ解空間上の距離が一定範囲外の第４の複数の解を一様抽出し、一様抽出した第４の複数の解に基づいて第４の最大評価値を再度推定する。これにより、探索された第１の最適解及び一定範囲外の解空間の第４の解候補からそれぞれ離れた解空間には、同等の最適解が存在しないことを確認することができる。
　本発明の更に他の態様に係る最適解判定方法において、解空間は、第１の制約条件における第１の解空間と第２の制約条件における第２の解空間とを含み、第１の解空間に属する第１の解候補を取得する第２４のステップと、第２の解空間上の複数の解を第５の複数の解として一様抽出する第２５のステップと、第５の複数の解のそれぞれに対応する第５の複数の評価値を取得する第２６のステップと、取得した第５の複数の評価値に基づいて、第５の複数の解の個数を超える個数の解を想定した場合の最大評価値を推定し、第５の最大評価値とする第２７のステップと、第２の解空間における解候補であって、第２４のステップで取得した第１の解候補からの解空間上の距離が近い解を近傍解として取得する第２８のステップと、近傍解の評価値を取得する第２９のステップと、近傍解の評価値が第５の最大評価値の信頼区間内に入るか否かを判定する第３０のステップと、を含むことが好ましい。
　本発明の更に他の態様によれば、第１の制約条件における第１の解空間における解候補の最適性を判定する場合、その解候補からの解空間上の距離が近い近傍解であって、第２の制約条件における第２の解空間における近傍解を取得し、近傍解の最適性を判定する。即ち、近傍解の最適性判定を代用し、近傍解が最適であれば、その近傍解に近い距離の第１の解空間における解候補も最適であると推定する。
　本発明の更に他の態様に係る最適解判定方法において、第１のステップは、組合せ最適化問題における組合せ可能なパターンのうち、組合せの一部により残りの組合せを考慮せずとも不適となることが確定するかどうかを識別することで、識別すべきパターンを削減するステップ、及び組合せ可能なパターンのうち、組合せの一部だけに差分があるパターン群の共通部分を抽出し、残りの組合せを共有することで、識別すべきパターンを削減するステップのうちの少なくとも一方を用いて、組合せ可能パターンを縮約して列挙索引化するデータ構造を用い、解空間上の解の総数を算出し、算出した総数以下の乱数を発生させ、発生させた乱数により特定されるパターンに対応する解を抽出することが好ましい。
　これにより、非常に大規模で、効率的な手段を用いてもなお全描出が難しい組合せ最適化問題であっても解を一様抽出することができ、解の最適性判定を行うことができる。
　本発明の更に他の態様に係る最適解判定方法において、組合せ最適化問題は、遺伝子制御ネットワークの組合せ最適化問題であることが好ましい。
　本発明の更に他の態様に係る最適解判定プログラムは、上記の最適解判定方法をコンピュータに実行させる。
　更に他の態様に係る発明は、組合せ最適化問題における解の最適性を判定する最適解判定装置であって、組合せ最適化問題の解空間上の複数の解を第１の複数の解として一様抽出する解抽出部と、一様抽出した第１の複数の解のそれぞれに対応する評価値を取得する第１の評価値取得部と、取得した第１の複数の評価値に基づいて、第１の複数の解の個数を超える個数の解を想定した場合の最大評価値を推定し、第１の最大評価値とする第１の最大評価値推定部と、解空間に属する解のうち少なくとも１つの解を解候補として取得する解取得部と、解候補に対応する評価値を第１の評価値取得部から取得し、取得した解候補に対応する評価値が第１の最大評価値の信頼区間内に入るか否かを判定し、解候補に対応する評価値が第１の最大評価値の信頼区間内に入ると判定された場合に、解候補を第１の最適解とする第１の判定部と、を備える。
　本発明の更に他の態様に係る最適解判定装置において、解空間上の解であって、第１の最適解からの解空間上の距離が一定範囲外の複数の解を第２の複数の解として一様抽出する第２の解抽出部と、第２の複数の解のそれぞれに対応する第２の複数の評価値を取得する第２の評価値取得部と、第２の複数の評価値に基づいて、第２の複数の解の個数を超える個数の解を想定した場合の最大評価値を推定し、第２の最大評価値とする第２の最大評価値推定部と、第１の最適解に対応する評価値を第２の評価値取得部から取得し、取得した第１の最適解に対応する評価値が、第２の最大評価値を超えているか否かを判定する第２の判定部と、を更に備えることが好ましい。
　本発明によれば、組合せ最適化問題の解の統計的な最適性判定が可能になり、組合せ最適化問題における解の最適性の判定を効率的かつ精度よく行うことができる。
ＲＮＡ発現行列データを示す図表及び遺伝子制御ネットワークを示す図 本発明の特徴を示す概念図 本発明に係る最適解判定装置のハードウェア構成を示すブロック図 本発明の第１の実施形態を示す機能ブロック図 要素「Ａ，Ｂ，Ｃ，Ｄ」からなる全体集合「Ｇ∧」を示す図 ４つの部分集合（Ｇ（１），Ｇ（２），Ｇ（３），Ｇ（４））を示す図 ３つの部分集合（Ｇ（１），Ｇ（３），Ｇ（４））の選択に対応する「パターン」を示す図 部分集合の選択に対応する「パターン」の３つの例と、各「パターン」が条件を満たしているか否かの判定結果とを示す図 図５から図８に示した集合分割問題における全てのパターン及び判定結果を、「２分グラフ」で網羅的に表現した図 フロンティア法の「枝刈り」により組合せ可能パターンが縮約される様子を示す図 フロンティア法の「節共有」により組合せ可能パターンが縮約される様子を示す図 フロンティア法の「枝刈り」及び「節共有」により組合せ可能パターンを縮約した結果を示す図 ノード＝｛Ａ，Ｂ，Ｃ｝で、エッジの数が３本のグラフ集合をＺＤＤ表現した図及び特定のパス（Ａ⇔Ｂ⇒Ｃ）に対応するグラフを示す図 ＺＤＤ表現したグラフ集合の全体のグラフ数（採用総数）の「数え上げ」を示す図 ＺＤＤ表現したグラフ集合から任意の指定番号のグラフを取り出す方法を示す図 本発明の第２の実施形態を示す機能ブロック図 ある遺伝子制御ネットワーク推定に対して本発明を適用した例を示すグラフ 本発明に係る最適解判定方法の第１の実施形態を示すフローチャート 本発明に係る最適解判定方法の第２の実施形態を示すフローチャートであり、特に最適十分性判定に成功した後に行われる最適必要性判定の処理を示すフローチャート 本発明に係る最適解判定方法の第４の実施形態の要部を示すフローチャート 本発明に係る最適解判定方法の第５の実施形態を示すフローチャート 制約条件Ｃにおける第１の解空間と制約条件Ｃ＋における第２の解空間とを模式的に表す図
　以下、添付図面に従って本発明に係る最適解判定方法、最適解判定プログラム及び最適解判定装置の好ましい実施の形態について説明する。
　＜本発明の概要＞
　組合せ最適化問題における解を探索する方法として、創薬分野に応用可能な遺伝子制御ネットワークを例に説明する。遺伝子制御ネットワークとは、遺伝子間の協調関係を有向グラフとして表現することで、例えば薬剤の作用機序などを読み解くための応用などが期待されている。
　まず、基本的な前提条件について説明する。
　(1) 図１は、ＲＮＡ発現行列データを示す図表及び遺伝子制御ネットワークを示す。
　図１において、Ａ、Ｂ、…、Ｚは遺伝子であり、Ｘ１、Ｘ２、…、Ｘｎはサンプルであり、遺伝子の数とサンプルの数との積だけデータが存在する。このＲＮＡ発現行列データを取得する。ＲＮＡ発現行列データは、ＮＧＳによるカバレッジデータであってもよいし、マイクロアレイによるシグナルデータでもよい。
　ＲＮＡ発現行列では、Ｍ個の細胞株等に対して、Ｎ個の遺伝子のＲＮＡ発現量が計測されており、データx(m,n)は、細胞株ｍにおける遺伝子ｎの発現量を示す。したがって、ＲＮＡ発現行列ＤはＭ×Ｎの数値行列データである。
　(2) 複数の遺伝子の関連性は、遺伝子制御ネットワークとして表すことができる。以下、遺伝子制御ネットワークをグラフＧとして表す。図１に示すようにグラフＧは、エッジ（矢印で示すノード（遺伝子）間の制御関係）の集合であり、例えばg1={(A,B),(A,C),(C,D)}は、「遺伝子AからB」「遺伝子AからC」「遺伝子CからD」の３つの制御関係が存在することを示す。グラフＧは、多数のサンプルでのＲＮＡ発現行列データＤから推定可能である。
　ただし、グラフＧには問題に応じて何らかの制約Ｃが課せられる。制約はモデル上、もしくは事前知識によって課せられる。例えばベイジアンネットワークモデルでは循環グラフを表現できないため、グラフＧは循環グラフであってはならない。（つまり、例えば「(A,B),(B,A)」や「(A,B),(B,C),(C,A)」といった部分集合を含んではいけない。）また、事前知識によってスケールフリーネットワーク性（ノードの次数分布がべき乗則に適合すること）が期待されていれば、そのような制約を設けることも考えられる。
　(3) 評価関数Ｓ（Ｄ，Ｇ）を用意する。これは、グラフＧがデータＤをどれだけ説明できているかを定量化したものである。例えば、前述のg1では、「遺伝子AからB」の制御関係に対しては、「x(m,B)=F(x(m,A))が当てはまるかどうか」を定量化する。Fは制御関係のモデル関数であり、定量化は、例えば罰則付き最大尤度（AIC（Akaike's Information Criterion）やBIC（Bayesian information criterion））が使われる。
　(4) 最適化問題を解く何らかのヒューリスティックな手法により、最適と思われるグラフＧ_1を獲得し、それに対する評価値Ｓ_1を取得する。
　上記は遺伝子制御ネットワークの推定の例であるが、(a) 何らかのデータＤに基づき、(b) 制約Ｃに従う集合Ｇを考え、(c) 評価関数Ｓの最大化（もしくは最小化）を試みて、(d) 特定の集合Ｇ_1を取得する、というところがポイントであって、かつ、組合せ最適化問題では一般的な構造である。
　したがって、本発明は、様々な他の組合せ最適化問題に対しても適用し得る。例えば、細胞の排他的かつ被覆的な遺伝子変異を探索する問題がある。これによって、例えば癌にとって重要な遺伝子変異又は作用機序を推定しようという応用が知られている。
　(1) Ｍ個の細胞に対して、Ｎ種類のＳＮＰ（Single Nucleotide Polymorphism）における変異データｘを取得する。データx(m,n)は、細胞mにおけるＳＮＰ変異nの有無を示す。ここで、ＳＮＰとはＤＮＡ（deoxyribonucleic acid）のうち変異が入り易い位置のことを言う。
　(2) 遺伝子座のセットを集合Ｇとして表す。集合ＧはＳＮＰの集合であり、例えばg1={1,3,4}は、「ＳＮＰ1,3,4」の３つのＳＮＰに注目することを示す。集合Ｇの要素はＮ種類のＳＮＰのいずれかであり、集合Ｇ全体でＳＮＰの組合せパターンに対応する。ただし、制約Ｃとして、集合ＧによってＭ個の細胞株を排他的に被覆することが要求される。即ち、「ＳＮＰ1に変異を持つ細胞の集合M1」「ＳＮＰ3に変異を持つ細胞の集合M3」「ＳＮＰ4に変異を持つ細胞の集合M4」を考えたとき、M1,M3,M4は互いに重複要素となる細胞を共有してはならず、かつ、M1,M3,M4全体ですべての細胞を網羅していなければならない。
　(3) 評価関数Ｓ（Ｇ）を用意する。これは、集合Ｇの何らかの特性を定量化したものである。例えば、ＳＮＰの集合Ｇを、予め用意した遺伝子制御ネットワークＦに載せたときの適合性などを評価することが考えられる。
　(4) 何らかのヒューリスティックな手法により、最適と思われるグラフＧ_1を獲得し、それに対する評価値Ｓ_1を取得する。
　上記(4)としては、例えば、グリーディ・ヒルクライミング法、焼きなまし法、タブーサーチ、遺伝的アルゴリズムなどが知られており、そのいずれを用いても良い。
　また、本発明はいわゆるバイオインフォマティクス以外の分野にも適用し得る。例えば遺伝子制御ネットワーク推定はベイジアンネットワークとして一般化されるので、多数の製品の様々な特性を測定してデータ化し、その特性同士の因果関係を推定する手法としても利用できる。組合せ最適化問題としては、例えばナップザック問題や巡回セールスマン問題などが知られ、様々な分野にて応用されており、本発明はそのいずれにも適用できる。
　さて、通常のヒューリスティック探索は、上記(4) でアルゴリズムは終了するので、ヒューリスティック探索により獲得したグラフＧ_1が、真の最適値かどうかを判断できなかった。
　例えば遺伝子制御ネットワーク問題の場合、グラフＧ_1がＲＮＡ発現行列データＤに対して真に最も適合しているのかどうかを判断できなかった。そのため、多額のコストを要する介入実験に踏み込むためには、例えばグラフＧ_1をバイオロジストが精査して妥当性を判断するなどの属人的で不確実な工程を要していた。
　そこで、本発明は、最適と思われるグラフＧ_1の最適性を見極めることができるようにする。
　図２は、本発明の特徴を示す概念図である。
　ヒューリスティック探索により探索した最適と思われるグラフＧ_1の最適性を見極めるために、ヒューリスティック探索により探索したグラフの抽出個数を超える個数を想定した場合の評価値（解空間全体の解の評価値（「スコア」とも言う））のうちの最大評価値（第１の最大評価値）Ｚを推定する。尚、第１の最大評価値Ｚの具体的な推定方法については後述する。
　続いて、最適と思われるグラフＧ_1（ローカル解）が、推定した第１の最大評価値Ｚの信頼区間内に入るか否かを判定する。そして、グラフＧ_1が第１の最大評価値Ｚの信頼区間内に入っていれば、グラフＧ_1は、探索空間（解空間）の全域での第１の最適解（グローバル解の一つ）であると判定できる。上記のような判定（解が十分条件を満たすか否かの最適十分性判定）が、本発明の特徴の一つである。
　ただし、上記の最適十分性判定の場合、ローカル解がグローバル解か否かの判定は可能であるが、唯一のグローバル解か否かの判定はできない。したがって、グラフＧ_1（ローカル解）の最適十分性判定に成功しても同等解が他にも存在する可能性があり、探索に「未練」が残る。
　本発明は、グラフＧ_1の最適十分性判定に成功した場合、そのグラフＧ_1からの解空間上の距離が一定範囲外の解空間（部分空間）上の解の評価値のうちの最大評価値（第２の最大評価値）Ｗを推定する。尚、第２の最大評価値Ｗの具体的な推定方法については後述する。
　続いて、グラフＧ_1が、推定した第２の最大評価値Ｗを超えているか否かの判定を更に行う。そして、グラフＧ_1が第２の最大評価値Ｗを超えていれば、部分空間には、グラフＧ_1と同等解が存在しないことになり、そのグラフＧ_1は唯一のグローバル解であると判定できる。上記のような判定（解が必要条件を満たすか否かの最適必要性判定）が、本発明の他の特徴の一つである。
　＜最適解判定装置＞
［装置構成］
　図３は本発明に係る最適解判定装置のハードウェア構成を示すブロック図である。
　図３に示す最適解判定装置１０は、コンピュータによって構成されており、主として各構成要素の動作を制御する中央処理装置（ＣＰＵ：Central Processing Unit）１２と、装置の制御プログラムが格納されたり、プログラム実行時の作業領域となる主メモリ１４と、液晶ディスプレイ、ＣＲＴディスプレイ等のモニタ装置２８の表示を制御するグラフィックボード１６と、ネットワーク５０と接続される通信インターフェース（通信Ｉ／Ｆ）１８と、本発明に係る最適解判定プログラムを含む各種のアプリケーションソフト、及び後述する最適解の判定結果等を保存するハードディスク装置２０と、ＣＤ－ＲＯＭドライブ２２と、キーボード３０のキー操作を検出して指示入力としてＣＰＵ１２に出力するキーボードコントローラ２４と、位置入力装置としてのマウス３２の状態を検出してモニタ装置２８上のマウスポインタの位置やマウス３２の状態等の信号をＣＰＵ１２に出力するマウスコントローラ２６とから構成されている。
　また、ネットワーク５０には、ＲＮＡ発現行列データを保存するデータベース４０が接続されている。ＲＮＡ発現行列データは、図１に示したように複数の細胞株（サンプル：Ｘ１、Ｘ２、…、Ｘｎ）における複数の遺伝子（Ａ、Ｂ、…、Ｚ）のＲＮＡ発現量を示す数値行列データである。また、ＲＮＡ発現量は、図示しないＮＧＳ（Next Generation Sequencer）などによりサンプルから取得されたものである。
　最適解判定装置１０は、通信インターフェース１８を介してデータベース４０にアクセスし、必要なＲＮＡ発現行列データを取得することができる。尚、ＲＮＡ発現行列データは、外部のデータベース４０に格納されたものを使用する場合に限らず、ＲＮＡ発現行列データをハードディスク装置２０に保存し、ハードディスク装置２０に保存されたＲＮＡ発現行列データを使用するようにしてもよい。
　［第１の実施形態］
　図４は、図３に示した最適解判定装置１０のＣＰＵ１２の機能を示す機能ブロック図であり、本発明の第１の実施形態を示す機能ブロック図である。
　ＣＰＵ１２は、ハードディスク装置２０に格納された最適解判定プログラムを実行することより各種の処理部として機能し、図４に示す第１の実施形態では、解抽出部１００、第１の評価値取得部及び第２の評価値取得部として機能する評価値取得部１０２、解入力部１０４、第１の最大評価値推定部１０６、第１の比較部１０８及び第１の判定部１１０としての機能を有する。
　解抽出部１００は、組合せ最適化問題（遺伝子制御ネットワーク）の解空間上の解（グラフ）を一様抽出する部分であり、本例では、ＺＤＤ（Zero-suppressed binary Decision Diagram）を用いたパス列挙索引化アルゴリズムによりグラフＧを列挙索引化する。制約Ｃの下でＺＤＤを構築することで、グラフＧの総数と、集合｛Ｇ｝の任意の要素を一様抽出できるようになる。遺伝子制御ネットワークでは、遺伝子をノード、制御関係をエッジとするグラフを考える。
　尚、ここでいう「解」は、「許容解（実行可能解）」を意味し、最適解とは限らないが、実行不可能ではない解のことを指す。すなわち（実行が不可能な）不適な解はあらかじめ排除されている。
　また、解抽出部１００は、「枝刈り」及び「節共有」の少なくとも一方を用いて、組合せ最適化問題における組合せ可能なパターンを縮約して列挙索引化するデータ構造とする。ここで、「枝刈り」とは、組合せ最適化問題における組合せ可能なパターンのうち、組合せの一部により残りの組合せを考慮せずとも不適となることが確定するかどうかを識別することで、識別すべきパターンを削減する処理をいう。また、フロンティア法の「節共有」とは、組合せ可能なパターンのうち、組合せの一部だけに差分があるパターン群の共通部分を抽出し、残りの組合せを共有することで、識別すべきパターンを削減する処理をいう。
　尚、これはＺＤＤだけに限定されるものではなく、組合せ最適化問題に応じて例えばＢＤＤ（Binary Decision Diagram）、あるいはπＤＤ（Permutation Decision Diagram）などのＺＤＤに類似する変形データ構造を用いても構わない。
　また、「枝刈り」の「不適かどうか」の判定には制約Ｃを用いる。例えば遺伝子制御ネットワークの場合では、採用済のエッジによって既に循環が発生すれば、残りのエッジを考慮しなくても不適であることが確定する。また、「節共有」の「共通かどうか」の判定においても、制約Ｃのもとでパターンのうち考慮済の部分と残りの要素とを勘案して、共通化できるかどうかを判断する。例えばエッジの数のみを考慮している場合、採用済のエッジの本数が同じであれば、「節共有」を適用できる。「枝刈り」及び「節共有」のアルゴリズムは、ＺＤＤの場合、フロンティア法として知られているので、それに従えばよい。
　尚、「枝刈り」のみを備えた手法として、分岐限定法によって解を列挙する手法を用いてもよい。ＺＤＤのうち「節共有」を無効化することで分岐限定法による列挙も可能となるし、「枝刈り」を無効化することも考えられる。しかし、望ましくはＺＤＤのように「枝刈り」及び「節共有」の両方とも利用することで、それによって効率的に解を列挙できる。
　解を一様抽出する手段としては、例えばランダム生成も考えられる。即ち、解候補をランダムに生成（例：グラフのエッジのありなしを乱数で決定）し、Ｇの制約を満たしていなければ再生成を繰り返す。制約条件が簡単であれば、ランダム生成の段階でその制約を織り込んでも構わない。例えばエッジ数を所定数以下にする場合、エッジありを選ぶ数に上限を設けておけばよい。一方、例えば循環グラフを禁止したい場合、単純なランダム生成において制約をかけるのは難しいので、ランダム生成されたグラフが循環しているかどうかを判定することが考えられる。しかし、理論的には、ランダム生成でも十分な実施数を重ねることで統計的推定は可能だが、特に判定及び再生成を繰り返す手法で、ランダム生成がカバーする解空間に対して制約下の解空間が小さい場合、例えば解空間のサイズを１：Ｎとすれば、１個の解を生成するのに平均してＮ個のランダム生成が必要になるので、効率が悪い。
　したがって、本手法では十分なサンプルサイズの確保が重要なため、特にＺＤＤの導入は効果が大きいと期待される。
　＜ＺＤＤの概説＞
　次に、ＺＤＤ及びフロンティア法について具体的に説明する。
　まず、組合せ最適化問題の一種である集合分割問題へのＺＤＤの適用を考える。
　集合分割問題は、ある全体集合に対する部分集合の列が与えられたとき、その幾つかを選んで、「選んだ部分集合同士に重複がない（相互排他）」、かつ「元の全体集合を尽くす（全体被覆）」のようなパターン（組合せ）を作れるか、という問題である。
　集合は、「要素の集まり」として定義される。図５に示すように「Ｇ∧」は全体集合を示し、「Ａ，Ｂ，Ｃ，Ｄ」が「要素」に対応する。
　要素を固定し、要素の有り又は無しを「１又は０」に割り振ると、例えば、図６に示すように部分集合が決まる。
　「部分集合を含むかどうか」を符号化すれば、図７に示すように「パターン」を表現できる。図７の例では、３つの部分集合（Ｇ（１），Ｇ（３），Ｇ（４））の選択に対応する「パターン」が示されている。
　続いて、パターンごとに「条件（被覆性及び排他性）を満たしているか」を判定する。
　図８には、部分集合の選択に対応する「パターン」の３つの例と、各「パターン」が条件を満たしているか否かの判定結果とが示されている。
　図８に示すように部分集合（Ｇ（１），Ｇ（３），Ｇ（４））の選択に対応する「パターン」は条件を満たし、部分集合（Ｇ（２），Ｇ（３），Ｇ（４））及び部分集合（Ｇ（３），Ｇ（４））の選択に対応する「パターン」は条件を満たさない。
　部分集合（Ｇ（２），Ｇ（３），Ｇ（４））は要素「Ｃ」が重複し、排他性を満たさず、部分集合（Ｇ（３），Ｇ（４））は要素「Ａ」が不足して被覆性を満たさないからである。
　さて、本例の集合分割問題における全てのパターン及び判定結果は、図９に示すように「２分グラフ」で網羅的に表現することができる。
　本例の部分集合は４個だったため、全てのパターンは、１６個（＝２４）であるが、部分集合がＮ個の場合、全てのパターンは、２Ｎとなる。「２分グラフ」の表現では、「２のべき乗」で枝及び葉が増えるという問題がある。
　組合せ問題は、「組合せ爆発」が発生し、有限時間で最適解を探索することができなくなることが知られているが、フロンティア法の「枝刈り」及び「節共有」の縮約技術により、所定の条件を満たす組合せ集合の効率的な（実用的な）「数え上げ」を実現することができる。
　図１０は、フロンティア法の「枝刈り」により組合せ可能パターンが縮約される様子を示す図である。
　図１０に示すように部分集合「Ｇ（１）」と「Ｇ（２）」とが同時に選択されると、要素「Ａ」が重複し、その時点で（後続選択に関わらず）不適が確定する。同様に部分集合「Ｇ（１）」と「Ｇ（２）」のどちらも選択しないと、要素「Ａ」を漏らすことになり、不適が確定する。そして、不適が確定すると、その時点で展開を打ち切り、「判定結果＝０」に即繋げる。
　このように「枝刈り」は、パターン選択の途中で不適が確定すると、その時点で展開を打ち切ることで組合せ可能パターンの縮約を図る。
　図１１は、フロンティア法の「節共有」により組合せ可能パターンが縮約される様子を示す図である。
　図１１に示すように２つの部分集合「Ｇ（１）」と「Ｇ（３）」とが選択される場合と、１つの部分集合「Ｇ（２）」が選択される場合とでは、いずれも「要素「Ａ」及び「Ｂ」を１回だけ含む」こと、かつ「要素「Ａ」及び「Ｂ」は、もはや含まないようにする」こと、が条件になるので、後続選択による採否判定が完全に一致する。
　そして、この場合には、別々に展開せずに、まとめて扱ってよい、つまり、同じ「節」を共有して構わない。
　このように「節共有」は、複数のパターン選択の後続が同一ならば、これらをまとめて扱うことで組合せ可能パターンの縮約を図る。
　図１２は、フロンティア法の「枝刈り」及び「節共有」により組合せ可能パターンを縮約した結果を示す図である。
　図１２に示すように、１６個のパターン（図９）に拡がるはずの枝及び葉を大幅に削減でき、網羅的に１つずつ判定する場合と完全に一致する判定結果を取得することができる。
　図１３は、ノード＝｛Ａ，Ｂ，Ｃ｝で、エッジの数が３本のグラフ集合をＺＤＤ表現した図である。
　図１３において、実線で示した「１枝」を通る場合のみ、そのエッジを含むグラフ、（点線で示した「０枝」もしくは飛ばしたエッジは含まない）であって、最終的に「１端」に到達したら、そのグラフを採用する。（「０端」に到達したグラフは採用しない。）
　ここで、例えば、（Ａ⇔Ｂ⇒Ｃ）のパスに対応するグラフは、図１３の左側の太線の矢印で示した経路で表される。
　図１４は、全体のグラフ数（採用総数）の「数え上げ」を示す図である。
　図１４に示すように採用総数は、判定結果を示す「１端」に「１」を付与し、「１端」から最上位のＺＤＤノードまで逆順に辿って数え上げることで算出することができる。
　数え上げは、下層ＺＤＤノードから各々の枝先の付与数を加算し、加算した数を自身に付与する。これを最上位のＺＤＤノードまで繰り返し、最上位のＺＤＤノードに付与された数値が、全体の採用総数（「１端」に到達するパスの総数）となる。本例の場合、全体の採用総数は、２０個になる。
　このように採用総数を算出する「数え上げ」は、ＺＤＤの重要な性質の一つである。
　図１５は、任意の指定番号のグラフを取り出す方法を示す図である。
　「数え上げ」後、任意の番号（本例の場合、１～２０の範囲内の番号）が指定されると、指定された番号（指定番号）にしたがって根から下ることで、指定番号に対応するグラフを取り出すことができる。
　例えば、「１２番」のグラフを取り出す場合、図１５の太線の矢印にしたがって最上位のノードから下る。まず、最上位のノード（Ａ，Ｂ）から「０枝」又は「１枝」のうち指定番号を含む枝に進む。本例では、「１枝」側の枝を進み、最上位のノード（Ａ，Ｂ）から下層のノード（Ａ，Ｃ）（図１４上で右側のノード（Ａ，Ｃ））に下る。「１枝」側の枝を進む場合、指定番号から「０枝」側の個数を引く。本例では、「１２番」のグラフを取り出すため、指定番号「１２」から「０枝」側の個数「１０」が引かれ、「２」になる。これを「１端」に到達するまで繰り返すことで、図１５の太線の矢印で示すパス（グラフ）を取り出すことができる。尚、図１５に示す「１２番」のグラフは、（Ａ⇒Ｂ⇔Ｃ）のパスに対応するグラフである。
　このように採用総数の「数え上げ」後、採用総数内の任意の番号を指定すると、指定した番号により一意に特定されるグラフを取り出すことができる。これにより採用総数以下の乱数を発生させることで、解空間上の解（グラフ）を「一様抽出」することができる。この解空間上の解の「一様抽出」は、ＺＤＤの重要な性質の一つである。
　図４に戻って、解抽出部１００は、ＺＤＤを用いたパス列挙索引化アルゴリズムにより遺伝子制御ネットワークの解空間上の解（グラフＧ）を一様抽出する。解空間上のグラフＧの総数は、ＺＤＤの重要な性質の一つである「数え上げ」により行うことができ、また、「数え上げ」により取得したグラフＧの総数以下の乱数（例えば、Ｍ系列（maximal length sequence）を用いた疑似乱数）を発生させ、各乱数により指定された指定番号に対応するグラフＧを取り出す（グラフＧを一様抽出する）。
　評価値取得部（第１の評価値取得部）１０２は、解抽出部１００により抽出されたグラフＧに評価値を付与する。例えば、グラフＧがＲＮＡ発現行列データＤをどれだけ説明できているかを定量化した評価関数Ｓ（Ｄ，Ｇ）を用意しておき、評価値取得部１０２は、抽出されたグラフＧに対応する評価値Ｓを評価関数Ｓ（Ｄ，Ｇ）に基づいて取得し、取得した評価値ＳをグラフＧに付与する。評価関数Ｓ（Ｄ，Ｇ）は、データベース４０に保存されたＲＮＡ発現行列データに基づいて評価値取得部１０２が作成してもよいし、予めＲＮＡ発現行列データに基づいて作成され、例えばデータベース４０に保存された評価関数Ｓ（Ｄ，Ｇ）を使用するようにしてもよい。
　解取得部１０４は、遺伝子制御ネットワークの最適化問題を解く何らかのヒューリスティックな手法により獲得された最適と思われる解（解候補）（以下、「グラフＧ_1」という）を取得する。また、入力したグラフＧ_1に対する評価値Ｓ_1は、評価値取得部１０２により取得することができる。尚、ヒューリスティックな手法としては、例えば、グリーディ法、ヒルクライミング法、焼きなまし法、タブーサーチ、遺伝的アルゴリズムなどが知られており、そのいずれを用いても良い。また、ヒューリスティックな探索は、本装置により行ってもよいし、外部の装置により行ってもよく、解取得部１０４は、何らかのヒューリスティックな手法により獲得された解候補（グラフＧ_1）を取得する。
　第１の最大評価値推定部１０６は、解空間上の解（グラフＧ）の最大評価値を推定する部分であり、本例では、解抽出部１００により一様抽出された複数の解（第１の複数の解）の評価値に基づいて、第１の複数の解の抽出個数を超える個数を想定した場合の評価値のうちの最大評価値（第１の最大評価値Ｚ）を推定する。
　具体的には、Ｕ，Ｖをそれぞれ自然数とすると、Ｕ×Ｖ個の解（第１の複数のグラフＧ）を一様抽出し、各々のグラフＧに評価値Ｓを与える。ここで、Ｕはブロックサイズであり、Ｖはブロック個数である。Ｕ、Ｖは、ある大きな数として設定する。例えばＵ、Ｖともに10,000に設定しても良く、この場合、一様抽出されるグラフＧの個数は、１億個（＝10,000×10,000)となる。
　第１の最大評価値推定部１０６は、Ｕ×Ｖ個のグラフＧをＶ個のブロックに分け、ブロック毎のＵ個のグラフＧの評価値のうちの区分最大値を取得する。したがって、区分最大値は、Ｖ個取得すことができる。そして、Ｖ個の区分最大値が一般極値分布（ＧＥＶ：generalized extreme value distribution）に従うものとして最大評価値（第１の最大評価値Ｚ）を最尤推定する。
　第１の最大評価値は、統計学的な裏付けを伴うものである。解空間上のグラフ｛Ｇ｝は本来は有限集合であり、厳密には縮退してしまうが、グラフＧの総数が十分に大きいため連続分布近似が適用できる。その場合、評価値Ｓに明らかに上限が存在するので、適切なＵ、Ｖの設定によってガンベル型になることが期待され、真の第１の最大評価値Ｚを信頼区間付で推定できる。
　第１の比較部１０８は、解取得部１０４が取得した解候補（グラフＧ_1）に対応する評価値Ｓ_1を評価値取得部１０２から取得し、取得した評価値Ｓ_1と推定した信頼区間付の第１の最大評価値Ｚとを比較する。
　第１の判定部１１０は、第１の比較部１０８による比較結果に基づいてグラフＧ_1の評価値Ｓ_1が、第１の最大評価値Ｚの信頼区間内に入るか否かを判定する。第１の最大評価値Ｚの信頼区間内に収まっていれば、グラフＧ_1は、解空間の全域での第１の最適解の一つ（「グラフＧ_1は十分」）であることが分かる。
　仮に、Ｚ≫Ｓ_1であれば、両者の評価値の差を解空間上の距離に変換し、現在推定しているグラフＧ_1が真の最適解（第１の最大評価値Ｚに対応する解）からどのくらい離れているかを推定してもよい。
　即ち、ヒューリスティック探索で探索された解候補（グラフＧ_1）が、真の最適値かどうかを判断できるようになる。
　［第２の実施形態］
　図１６は、図３に示した最適解判定装置１０のＣＰＵ１２の機能を示す機能ブロック図であり、本発明の第２の実施形態を示す機能ブロック図である。尚、図１６において、図４に示した第１の実施形態と共通する部分には同一の符号を付し、その詳細な説明は省略する。
　図１６に示す第２の実施形態はて第１の実施形態と比較して第２の最大評価値推定部１１２、第２の比較部１１４及び第２の判定部１１６が主として追加されている。
　第１の実施形態では、解候補（Ｇ_1）の最適性十分性判定に成功しても、解候補（Ｇ_1）以外にも同等に確からしい別の解候補が存在する可能性は排除できない。
　第２の実施形態は、解候補（Ｇ_1）以外には同等に確からしい別の解候補が存在しないことを判定可能にするものである。
　図１６において、第１の実施形態により解候補（グラフＧ_1）の最適性十分性判定に成功すると、その後、解抽出部１００は、グラフＧ_1からの解空間上の距離が一定範囲外の解を列挙索引化するＺＤＤを構築する。このＺＤＤの構築は、ＺＤＤを構築する際のフロンティア法において、構築途上のグラフＧとグラフＧ_1との共通又は非共通エッジをカウントすることで実現できる。解抽出部１００は、再構築したＺＤＤを用いたパス列挙索引化アルゴリズムにより遺伝子制御ネットワークの解空間上の第２の複数の解（グラフＧ）であって、グラフＧ_1から解空間上の距離が一定範囲外の第２の複数のグラフＧを一様抽出する。
　第２の最大評価値推定部１１２は、一様抽出された第２の複数のグラフにそれぞれ対応する評価値に基づいて、第２の複数の解の抽出個数を超える個数を想定した場合の評価値のうちの最大評価値（第２の最大評価値Ｗ）を推定する。
　具体的には、第１の最大評価値推定部１０６による第１の最大評価値Ｚの推定と同様の手法により第２の最大評価値Ｗを推定する。即ち、Ｐ，Ｑをそれぞれ自然数とすると、Ｐ×Ｑ個の解（第２の複数のグラフＧ）を一様抽出し、各々のグラフＧに評価値Ｓを与える。ここで、Ｐはブロックサイズであり、Ｑはブロック個数である。Ｐ、Ｑは、第１の最大評価値Ｚを推定する際に一様抽出したＵ、Ｖと同じでもよいし、異なっていてもよい。
　第２の最大評価値推定部１１２は、Ｐ×Ｑ個のグラフＧをＱ個のブロックに分け、ブロック毎のＵ個のグラフＧの評価値のうちの区分最大値を取得する。したがって、区分最大値は、Ｑ個取得すことができる。そして、Ｑ個の区分最大値が一般極値分布に従うものとして第２の最大評価値Ｗを最尤推定する。
　第２の比較部１１４は、解取得部１０４が取得した解候補（グラフＧ_1）に対応する評価値Ｓ_1を評価値取得部（第２の評価値取得部）１０２から取得し、取得した評価値Ｓ_1と推定した信頼区間付の第２の最大評価値Ｗとを比較する。
　第２の判定部１１６は、第２の比較部１１４による比較結果に基づいてグラフＧ_1の評価値Ｓ_1が、第２の最大評価値Ｗ（信頼区間付の第２の最大評価値Ｗの範囲）を超えているか否かを判定する。
　第１の最大評価値Ｚは、解空間全体の最大値を推定したものであり、第２の最大評価値Ｗは、部分空間の最大値を推定したものなので、基本的にはＷ≦Ｚは明らかである（サンプルサイズ等によっては確率的にＷ＞Ｚとなる場合もある）。
　その上で、Ｓ_1≫Ｗ（信頼区間を外れる等）であれば、グラフＧ_1は第１の最大評価値を与えるグラフであって、しかもグラフＧ_1から解空間上で離れた範囲には、同等以上にＲＮＡ発現行列データＤを説明できるグラフ構造はないと判断できる。
　尚、グラフＧ_1から離間させる一定範囲の距離は、事前に設定する必要があるが、これはグラフの特性から設定してもよいし、経験的に設定してもよいし、例えばＳ_1≫Ｗになるところまで徐々に離れる距離を大きくしていってもよい。例えば、十分大きな距離の設定値から２分検索などによって効率的に適切な距離を繰り返し探索しても構わない。
　言うまでもないことだが、距離ゼロであれば、Ｗ≒Ｚとなり、第１の実施形態と変わらない結果になろうし、ゼロ以外の最短の距離の場合はグラフＧ_1のみしか排除しないため、ある程度の距離を設定しないと、結果には大きな違いは出にくいと想定される。
　これにより、第２の最大評価値Ｗを超えていれば、グラフＧ_1は、解空間の全域での唯一の最適解（「グラフＧ_1は必要十分」）であることが分かる。即ち、ヒューリスティック探索で得られたグラフＧ_1以外に、最適値に相当するグラフが存在しないと判断できるようになる。
　尚、第２の実施形態は、第１の実施形態において、グラフＧ_1の最適十分性の判定に成功した場合に、続けて行うのが通例であるが、何らかの理由によって決められたグラフＧが存在する場合、そのグラフＧから離れた範囲の解候補とグラフＧとを直接比較する手段として利用されてもよい。
　従来のヒューリスティック探索では、例えば初期値をランダムに変えるとか、データにノイズを与えるなどの工夫で繰り返し探索を行うなどの方法もあったが、それ自体もヒューリスティックな判定法であるのに対して、本発明は統計的根拠に裏付けられた手法である。
　図１７は、ある遺伝子制御ネットワーク推定に対して本発明を適用した例を示すグラフである。
　図１７に示すグラフの横軸は乖離度（距離）であり、縦軸は評価値である。点線は、獲得グラフに対する到達評価値（第１の最大評価値Ｚ）である。また、白抜きの線は、各乖離度に対する最適値の推定値である。
　尚、ここでいう「乖離度」とは一般に距離またはそれに準じる指標である。例えば、集合をバイナリ配列として表現する場合はハミング距離としてもよい。問題に応じて、編集距離などの他の距離指標、あるいは、準距離、半距離等の指標を用いてもよい。また、それらに近しい乖離度指標を定義してもよい。
　乖離度ゼロでの最適値推定範囲は、到達評価値を包含しているため、到達評価値は最適値であると判定された。
　乖離度を増やすに連れて推定範囲に対応する値は徐々に低下し、乖離度５～６で推定範囲（第２の最大評価値Ｗ）が到達評価値を外れたので、乖離度５～６以上の範囲には獲得グラフ同等以上の評価値を有するグラフは存在しないことと判定された。
　また、実線はヒューリスティック探索で得られた解候補に対応する実測値を示すが、本発明が実態を正しく推定できていることを示す。
　＜最適解判定方法＞
　［第１の実施形態］
　図１８は、本発明に係る最適解判定方法の第１の実施形態を示すフローチャートである。
　図１８において、図４に示した解抽出部１００は、組合せ最適化問題の一つある遺伝子制御ネットワークの解空間上の解（グラフＧ）を一様抽出する（ステップＳ１０（第１のステップ））。本例では、前述したようにＺＤＤを用いたパス列挙索引化アルゴリズムによりグラフＧを一様抽出する。
　続いて、評価値取得部１０２は、一様抽出されたグラフＧに対して評価値Ｓを付与する（ステップＳ１２）。例えば、グラフＧがＲＮＡ発現行列データＤをどれだけ説明できているかを定量化した評価関数Ｓ（Ｄ，Ｇ）を用意しておき、評価値取得部１０２は、一様抽出されたグラフＧに対応する評価値Ｓを評価関数Ｓ（Ｄ，Ｇ）に基づいて取得し、取得した評価値ＳをグラフＧに付与する。
　第１の最大評価値推定部１０６は、一様抽出された複数のグラフＧ（第１の複数の解）に対して、ステップＳ１２により付与された第１の複数の評価値Ｓを取得し（ステップＳ１４、第２のステップ）、取得した第１の複数の評価値Ｓに基づいて、第１の複数のグラフＧの抽出個数を超える個数を想定した場合の評価値のうちの最大評価値（第１の最大評価値）Ｚを推定する（ステップＳ１６、第３のステップ）。
　具体的には、第１の複数のグラフＧは、Ｕ，Ｖをそれぞれ自然数とすると、Ｕ×Ｖ個のグラフであり、第１の最大評価値推定部１０６は、Ｕ×Ｖ個のグラフＧをＶ個のブロックに分け、ブロック毎のＵ個のグラフＧの評価値の区分最大値をＶ個取得し、Ｖ個の区分最大値を用いて、区分最大値が一般極値分布に従うものとして最大評価値（信頼区分付の第１の最大評価値Ｚ）を最尤推定する。
　続いて、解取得部１０４は、遺伝子制御ネットワークの最適化問題を解く何らかのヒューリスティックな手法により獲得された最適と思われる解候補（グラフＧ_1）を取得する（ステップＳ１８、第４のステップ）。ヒューリスティックな手法としては、例えば、グリーディ法、ヒルクライミング法、焼きなまし法、タブーサーチ、遺伝的アルゴリズムなどが知られており、そのいずれを用いても良い。
　評価値取得部１０２（第１の評価値取得部）により、取得したグラフＧ_1に対する評価値Ｓ_1を取得する（ステップＳ２０、第５のステップ）。
　第１の比較部１０８は、ステップＳ２０で取得した評価値Ｓ_1と、ステップＳ１６で推定した信頼区間付の第１の最大評価値Ｚとを比較する（ステップＳ２２、第６のステップ）。
　第１の判定部１１０は、第１の比較部１０８による比較結果に基づいてグラフＧ_1の評価値Ｓ_1が、信頼区間付の第１の最大評価値Ｚの範囲に入るか否かを判定する（ステップＳ２４、第７のステップ）。即ち、第１の判定部１１０は、評価値Ｓ_1が信頼区間付の第１の最大評価値の信頼区間内に入っている場合には、グラフＧ_1は最適解としての十分条件を満たしている（最適十分性あり）と判定する。
　第１の判定部１１０による判定結果は、図３に示したモニタ装置２８に表示され、又はハードディスク装置２０に保存され、又は図示しないプリンタにプリント出力される（ステップＳ２６、第８のステップ）。尚、第１の判定部１１０による判定結果は、最適十分性の有無に限らず、最適十分性がない場合には、評価値Ｓ_1と第１の最大評価値Ｚとの差を解空間上の距離に変換し、現在推定しているグラフＧ_1が真の最適解（第１の最大評価値Ｚに対応する解）からどのくらい離れているかを判定してもよい。
　第１の実施形態によれば、ヒューリスティック探索で探索された解候補（グラフＧ_1）が、真の最適値としての最適十分性を有するか否かを判定することができる。また、最適十分性の判定に失敗した場合でも、ヒューリスティック探索が不十分なために、他の最適解が存在することが示唆されたことになり、更に失敗程度によって、ある程度は最適解に近いことを主張したり、あるいは、最適解は別にあることを留意した上で、解候補を利用したりしても良い。即ち、最適十分性判定の成否に関わらず、最適十分性判定の情報は有用である。
　［第1の実施形態の変形例］
　本手順の最適十分性判定に失敗した場合、ヒューリスティックな同じ探索を異なる設定などで繰り返しても良いが、予めヒューリスティック探索を行う複数の探索法（第１の探索法、第２の探索法等）を用意し、複数の探索法を切り替えて使用する方法が考えられる。
　例えば、ヒューリスティックな探索法として、演算コストは小さい（探索時間は短い）が解の精度が低い第１の探索法と、第１の探索法よりも演算コストは大きい（探索時間は長い）が解の精度が高い第２の探索法とを準備しておき、最初に第１の探索法により探索した第１の解候補（グラフＧ_1）の最適十分性を判定し、最適十分性の判定に失敗した場合のみ、第２の探索法により探索した第２の解候補（グラフＧ_1）の最適十分性を判定する。
　これらの探索法の切り替えは、ヒューリスティックな探索法同士で切り替えてもよいし、近似性の保証がある程度ある近似アルゴリズムや厳密解を求める手法などへ切り替えてもよい。また、探索法は３つ以上用意して順次切り替えても構わない。探索法の切り替えは、方法自体に限らず、同一の探索法の収束判定等によって実現しても構わない。例えば、繰り返しサーチによって精度を高める探索法において、所定回の探索結果を第１の最大評価値Ｚで判定し、第１の最大評価値Ｚの信頼区間内に達するまで探索を繰り返しても良い。
　また、この場合、先行して最適十分性を判定するための第１の最大評価値Ｚを取得しておいても構わない。
　［第２の実施形態］
　図１９は、本発明に係る最適解判定方法の第２の実施形態を示すフローチャートであり、特に図１８に示した第１の実施形態による最適十分性判定に成功した後に行われる最適必要性判定の処理に関して示している。
　図１９において、図１６に示した解抽出部１００（第２の解抽出部）は、遺伝子制御ネットワークの解空間上の第２の複数のグラフＧであって、最適十分性判定に成功したグラフＧ_1からの解空間上の距離が一定範囲外の複数（第２の複数）のグラフＧを一様抽出する（ステップＳ３０、第９のステップ）。尚、第２の複数のグラフＧの一様抽出は、グラフＧ_1からの解空間上の距離が一定範囲外の解を列挙索引化するＺＤＤを構築し、構築したＺＤＤを用いたパス列挙索引化アルゴリズムにより行うことができる。
　続いて、評価値取得部１０２は、一様抽出された第２の複数のグラフＧに対してそれぞれ評価値Ｓを取得する（ステップＳ３２、第１０のステップ）。グラフＧに対する評価値Ｓの取得は、図１８に示した第１の実施形態のステップＳ１２と同様に行うことができる。
　次に、第２の最大評価値推定部１１２は、ステップＳ３２により取得した第２の複数の評価値Ｓに基づいて、第２の複数のグラフＧの抽出個数を超える個数を想定した場合の評価値のうちの最大評価値（第２の最大評価値）Ｗを推定する（ステップＳ３４、第１１のステップ）。
　具体的には、第２の複数のグラフＧは、Ｐ，Ｑをそれぞれ自然数とすると、Ｐ×Ｑ個のグラフであり、第２の最大評価値推定部１１２は、Ｐ×Ｑ個のグラフＧをＱ個のブロックに分け、ブロック毎のＰ個のグラフＧの評価値の区分最大値をＱ個取得し、Ｑ個の区分最大値を用いて、区分最大値が一般極値分布に従うものとして最大評価値（信頼区分付の第２の最大評価値Ｗ）を最尤推定する。
　続いて、第２の比較部１１４は、最適十分性判定に成功したグラフＧ_1の評価値Ｓ_1と、ステップＳ３４で推定した信頼区間付の第２の最大評価値Ｗとを比較する（ステップＳ３６、第１２のステップ）。
　第２の判定部１１６は、第２の比較部１１４による比較結果に基づいてグラフＧ_1の評価値Ｓ_1が、信頼区間付の第２の最大評価値Ｗの範囲を超えているか否かを判定する（ステップＳ３８、第１３のステップ）。即ち、第２の判定部１１６は、評価値Ｓ_1が信頼区間付の第２の最大評価値の範囲を超えている場合には、グラフＧ_1と同等解は他には存在せず、グラフＧ_1は唯一の最適解としての必要条件を満たしている（最適必要性あり）と判定する。
　第２の判定部１１６による判定結果は、図３に示したモニタ装置２８に表示され、又はハードディスク装置２０に保存され、又は図示しないプリンタにプリント出力される（ステップＳ４０、第１４のステップ）。
　第２の実施形態によれば、ヒューリスティック探索で探索された解候補（グラフＧ_1）が、真の最適値としての最適十分性判定に成功すると、更にグラフＧ_1の最適必要性を判定するため、探索されたグラフＧ_1から離れた解空間には、同等の評価値を有する解が存在しないことを確認することができる。
　［第３の実施形態］
　本発明に係る最適解判定方法の第３の実施形態は、図１９に示した第２の実施形態において、最適十分性判定に失敗した場合の処理を含む。
　即ち、最適十分性判定に失敗した場合、図１９に示したステップＳ３０（第１５のステップ）において、グラフＧ_1からの一定範囲を解空間上で拡大し、拡大した一定範囲外の解空間上で第３の複数のグラフＧ_3（第３の複数の解）を一様抽出する。
　そして、拡大した一定範囲外の解空間上で抽出した第３の複数のグラフＧに対応する第３の複数の評価値の取得（ステップＳ３２、第１６のステップ）、第３の複数の評価値に基づいて第３の複数のグラフＧの個数を超える個数の解を想定した場合の第３の最大評価値Ｗの推定（ステップＳ３４、第１７のステップ）、グラフＧ_1の評価値Ｓ_1と第３の最大評価値Ｗとの比較（ステップＳ３６）、及びグラフＧ_1の最適必要性の判定（ステップＳ３８、第１８のステップ）等を再度実行する。
　最適十分性判定に失敗した場合には、一定範囲を拡大させる上記の処理を最適十分性判定に成功するまで、一定範囲を徐々に拡大して複数回繰り返してもよい。
　［第４の実施形態］
　本発明に係る最適解判定方法の第４の実施形態は、図１９に示した第２の実施形態において、最適十分性判定に失敗した場合の他の処理を含む。
　最適十分性判定に失敗した場合（即ち、グラフＧ_1の評価値Ｓ_1が第２の最大評価値Ｗを超えていないとの判定結果が出力された場合）、図２０に示すようにグラフＧ_1からの解空間上の距離が一定距離離れた第４の解候補（グラフＧ_2)であって、第２の最大評価値Ｗに対応する評価値Ｓ_2を有するグラフＧ_2を取得する(ステップＳ４２、第１９のステップ)。
　グラフＧ_2の取得は、グラフＧ_1からの距離が一定範囲外の解空間の範囲で解候補（グラフＧ_2)のヒューリスティック探索を実行し、グラフＧ_2の評価値Ｓ_2が、Ｗ≒Ｚと同等となるグラフＧ_2を取得することにより行う。
　次に、図１９に示したステップＳ３０の代わりに、グラフＧ_1及びグラフＧ_2のそれぞれからの解空間上の距離が一定範囲外の第４の複数の解を一様抽出する（ステップＳ４４、第２０のステップ）。
　その後、図１９に示したステップＳ３２に遷移させ、上記のようにして拡大した一定範囲外の解空間上で抽出した第４の複数のグラフＧに対応する第４の複数の評価値の取得（ステップＳ３２、第２１のステップ）、第４の複数の評価値に基づいて第４の複数のグラフＧの個数を超える個数の解を想定した場合の第４の最大評価値Ｗの推定（ステップＳ３４、第２２のステップ）、グラフＧ_1の評価値Ｓ_1と第２の最大評価値Ｗとの比較（ステップＳ３６）、及びグラフＧ_1の最適必要性の判定（ステップＳ３８、第２３のステップ）等を再度実行する。
　そして、最適十分性判定に再度失敗した場合には、新たにグラフＧ_3等を追加し、新たな第４の最大評価値Ｗを推定して最適必要性の判定を繰り返し行う。
　尚、一定範囲（意味のある距離）は、組合せ最適化問題のもとの要請から決まる。例えば遺伝子制御ネットワークの場合、距離ｄは異なるエッジの本数を意味するので、作用機序解明等で許容されるエッジの間違い数によって解釈すればよい。エッジの総数がＮ本程度と見込まれる場合、ｄ／Ｎで誤答率を表す指標も考えられるので、例えば誤答率５％を許容してＮ=100が見込まれる場合は、ｄ=5等と設定すればよい。
　［第５の実施形態］
　解空間全体から解を一様抽出することができない場合がある。解空間全体のＺＤＤが構築できない事例等、そもそも一様抽出手段が確保できない場合である。
　その場合に適用可能な本発明に係る最適解判定方法の第５の実施形態について説明する。
　図２１は本発明に係る最適解判定方法の第５の実施形態を示すフローチャートである。
　まず、制約条件ＣだけではＺＤＤを構築できない場合、ＺＤＤを構築できる、より厳しい制約条件Ｃ＋を考える。例えばグラフ問題では、制約条件Ｃが「非循環」のみであれば、「全域林」を付与すること等が考えられる。組合せ最適化問題の解空間は、上記の制約条件Ｃ（第１の制約条件）における第１の解空間と、制約条件Ｃ＋（第２の制約条件）における第２の解空間とを含むものとする。
　図２１において、第１の解空間に属する解(制約条件Ｃにおける｛Ｇ｝）の中から第１の解候補（グラフＧ_1）を取得する（ステップＳ５０、第２４のステップ）。第１の解空間では、ＺＤＤ等を構築できていなくても、ヒューリスティックな手法等でグラフＧ_1の探索は可能である。
　続いて、第２の解空間上の解（制約条件Ｃ＋におけるグラフ｛Ｇ＋｝）を列挙し、ステップＳ５０で取得したグラフＧ_1の、制約条件Ｃ＋における近傍解Ｇ_1＋を探索する（ステップＳ５２、第２８のステップ）。例えば、制約条件Ｃ＋におけるＺＤＤは記述できることから、制約条件Ｃ＋に対してさらに「Ｇ_1からの解空間上の距離が一定範囲内であること」を付け加えてＺＤＤを構築し、そのうち最小距離のものを選ぶことなどにより近傍解Ｇ_1＋を探索できる。
　次に、近傍解Ｇ_1＋の評価値Ｓ_1＋を取得する（ステップＳ５４、第２９のステップ）。近傍解Ｇ_1＋に対する評価値Ｓ_1＋の取得は、図１８に示した第１の実施形態のステップＳ１２と同様に行うことができる。
　ステップＳ５４で取得した近傍解Ｇ_1＋の評価値Ｓ_1＋と、信頼区間付の第５の最大評価値Ｚとを比較する（ステップＳ５６）。尚、第５の最大評価値Ｚは、第２の解空間の制約条件Ｃ＋におけるグラフ｛Ｇ＋｝の一様抽出に基づいて前述した手法によって推定することできる。即ち、第２の解空間上の複数の解を第５の複数の解（制約条件Ｃ＋におけるグラフ｛Ｇ＋｝）として一様抽出し（第２５のステップ）、第５の複数の解のそれぞれに対応する第５の複数の評価値を取得し（第２６のステップ）、取得した第５の複数の評価値に基づいて、第５の複数の解の個数を超える個数の解を想定した場合の第５の最大評価値の推定する（第２７のステップ）。
　そして、近傍解Ｇ_1＋の評価値Ｓ_1＋と、信頼区間付の第１の最大評価値Ｚとの比較結果に基づいて、グラフＧ_1の最適十分性を判定する（ステップＳ５８、第３０のステップ）。即ち、評価値Ｓ_1＋が第１の最大評価値Ｚの信頼区間内に入るか否かを判定し、評価値Ｓ_1＋が第１の最大評価値Ｚの信頼区間内に入り、近傍解Ｇ_1＋の最適十分性があると判定されれば、グラフＧ_1も最適十分性があると判定する。近傍解Ｇ_1＋が最適なら、少なくとも「より制約が厳しい条件を課した場合の最適解から、最も近いところに解を見つけている」ことは担保できるからである。
　図２２は、制約条件Ｃにおける第１の解空間と、制約条件Ｃ＋における第２の解空間とを模式的に表す図である。ＺＤＤを構築できない制約条件Ｃにおける第１の解空間上でヒューリスティックな手法等で探索したグラフＧ_1の最適性判定を行う場合、ＺＤＤを構築できる、より厳しい制約条件Ｃ＋における第２の解空間上で近傍解Ｇ_1＋（グラフＧ_1からの解空間上の距離が一定範囲内の最小距離のもの）を探索し、この近傍解Ｇ_1＋の最適性判定を代用し、グラフＧ_1の最適性判定を行う。
　第２の解空間は、制約が追加されるので「Ｃ＋⊆Ｃ」となる（図２２）。なぜなら、制約条件Ｃに含まれる制約は、モデルの前提条件であったり、事前知識による妥当な制約であったりするため、基本的には外すべきではないからである。言うまでもなく、追加される制約は、事前知識等によって問題に対してある程度の妥当性が想定されることが望ましい（が、必ずしも完全な妥当性が担保されていなくても止むを得ない）。
　しかし、もし制約を外した方が探索容易で、かつ、制約を外しても妥当性の致命的な毀損には至らない場合であれば、「Ｃ⊆Ｃ＋」になるような制約、即ち制約の解除を考えても構わない。
　いずれにせよ、制約条件ＣとＣ＋との差及びグラフＧ_1と近傍解Ｇ_1＋との最小距離ｄは、代用した近傍解Ｇ_1＋と実際に利用する到達解との乖離度を示す。近傍解Ｇ_1＋の探索については、例えばヒューリスティックな手法によって、到達解の最小距離ｄの近傍については探索していると言えるのであれば、実用上は問題にならないと考えられる。一方、制約条件ＣとＣ＋との違いの影響は、それによって評価値がどの程度大きく変化するかに依存する。評価値にそれほど大きな影響を与えないような制約条件で、最小距離dが妥当な範囲であれば、近傍解Ｇ_1＋の最適性判定を代用し、グラフＧ_1の最適性判定を行う信頼度は高いと考えてよい。以上から、より広い組合せ最適化問題に対して本発明の効果を及ぼすことができるようになる。
　［その他］
　本実施形態の最適解判定装置１０は、例示に過ぎず、他の構成に対しても本発明を適用することが可能である。各機能構成は、任意のハードウェア、ソフトウェア、或いは両者の組合せによって適宜実現可能である。例えば、上述の最適解判定装置１０の各部における処理をコンピュータに実行させる最適解判定プログラム、そのような最適解判定プログラムを記録したコンピュータ読み取り可能な記録媒体（非一時的記録媒体）に対しても、本発明を適用することが可能である。
　また、本実施形態において、例えば、解抽出部１００、評価値取得部１０２、解取得部１０４、第１の最大評価値推定部１０６、第１の比較部１０８、及び第１の判定部１１０等の各種の処理を実行する処理部（processing unit）のハードウェア的な構造は、次に示すような各種のプロセッサ（processor）である。各種のプロセッサには、ソフトウェア（プログラム）を実行して各種の処理部として機能する汎用的なプロセッサであるＣＰＵ（Central Processing Unit）、ＦＰＧＡ（Field Programmable Gate Array）などの製造後に回路構成を変更可能なプロセッサであるプログラマブルロジックデバイス（Programmable Logic Device：ＰＬＤ）、ＡＳＩＣ（Application Specific Integrated Circuit）などの特定の処理を実行させるために専用に設計された回路構成を有するプロセッサである専用電気回路などが含まれる。
　１つの処理部は、これら各種のプロセッサのうちの１つで構成されていてもよいし、同種または異種の２つ以上のプロセッサ（例えば、複数のＦＰＧＡ、あるいはＣＰＵとＦＰＧＡの組み合わせ）で構成されてもよい。また、複数の処理部を１つのプロセッサで構成してもよい。複数の処理部を１つのプロセッサで構成する例としては、第１に、クライアントやサーバなどのコンピュータに代表されるように、１つ以上のＣＰＵとソフトウェアの組合せで１つのプロセッサを構成し、このプロセッサが複数の処理部として機能する形態がある。第２に、システムオンチップ（System On Chip：ＳｏＣ）などに代表されるように、複数の処理部を含むシステム全体の機能を１つのＩＣ（Integrated Circuit）チップで実現するプロセッサを使用する形態がある。このように、各種の処理部は、ハードウェア的な構造として、上記各種のプロセッサを１つ以上用いて構成される。
　さらに、これらの各種のプロセッサのハードウェア的な構造は、より具体的には、半導体素子などの回路素子を組み合わせた電気回路（circuitry）である。
　また、本発明は、プロセッサを有する最適解判定装置であって、プロセッサが、組合せ最適化問題の解空間上の複数の解を第１の複数の解として一様抽出し、一様抽出した第１の複数の解のそれぞれに対応する第１の複数の評価値を取得し、取得した第１の複数の評価値に基づいて、第１の複数の解の個数を超える個数の解を想定した場合の最大評価値を第１の最大評価値として推定し、解空間に属する解のうち少なくとも１つの解を解候補として取得し、取得した解候補に対応する評価値が第１の最大評価値の信頼区間内に入るか否かを判定する最適解判定装置を含む。
　更に、本発明は上述した実施形態に限定されず、本発明の精神を逸脱しない範囲で種々の変形が可能であることは言うまでもない。
１０　最適解判定装置
１２　ＣＰＵ
１４　主メモリ
１６　グラフィックボード
１８　通信インターフェース
２０　ハードディスク装置
２２　ＣＤ－ＲＯＭドライブ
２４　キーボードコントローラ
２６　マウスコントローラ
２８　モニタ装置
３０　キーボード
３２　マウス
４０　データベース
５０　ネットワーク
１００　解抽出部
１０２　評価値取得部
１０４　解入力部
１０６　第１の最大評価値推定部
１０８　第１の比較部
１１０　第１の判定部
１１２　第２の最大評価値推定部
１１４　第２の比較部
１１６　第２の判定部
Ｚ　第１の最大評価値
Ｗ　第２の最大評価値

Claims (16)

  1. 　組合せ最適化問題における解の最適性をコンピュータにより判定する最適解判定方法であって、
    　前記組合せ最適化問題の解空間上の複数の解を第１の複数の解として一様抽出する第１のステップと、
    　前記第１のステップにより一様抽出した前記第１の複数の解のそれぞれに対応する第１の複数の評価値を取得する第２のステップと、
    　取得した前記第１の複数の評価値に基づいて、前記第１の複数の解の個数を超える個数の解を想定した場合の最大評価値を推定し、第１の最大評価値とする第３のステップと、
    　前記解空間に属する解のうち少なくとも１つの解を解候補として取得する第４のステップと、
    　前記解候補に対応する評価値を取得する第５のステップと、
    　前記第５のステップにより取得した前記解候補に対応する評価値が前記第１の最大評価値の信頼区間内に入るか否かを判定する第６のステップと、
    　前記第６のステップにおいて前記解候補に対応する評価値が前記第１の最大評価値の信頼区間内に入ると判定された場合に、前記解候補を第１の最適解とする第７のステップと、
    　を含む最適解判定方法。
  2. 　前記第１の複数の解は、Ｕ，Ｖをそれぞれ自然数とすると、Ｕ×Ｖ個の解であり、
    　前記第３のステップは、前記Ｕ×Ｖ個の解をＶ個のブロックに分け、前記ブロック毎にＵ個の解の評価値の区分最大値をＶ個取得し、前記Ｖ個の区分最大値を用いて、前記区分最大値が一般極値分布に従うものとして前記第１の最大評価値を推定する、
    　請求項１に記載の最適解判定方法。
  3. 　前記第７のステップによる判定結果を出力する第８のステップを更に含む、
    　請求項１又は２に記載の最適解判定方法。
  4. 　演算コストは小さいが解の精度が低い第１の探索法と、前記第１の探索法よりも演算コストは大きいが解の精度が高い第２の探索法とを有し、
    　前記第４のステップは、最初に前記第１の探索法により探索された第１の解候補を入力し、前記第１の解候補の評価値が前記第１の最大評価値の範囲に入らない場合のみ、前記第２の探索法により探索された第２の解候補を入力する、
    　請求項１から３のいずれか１項に記載の最適解判定方法。
  5. 　前記解空間上の解であって、前記第１の最適解からの前記解空間上の距離が一定範囲外の複数の解を第２の複数の解として一様抽出する第９のステップと、
    　前記第２の複数の解のそれぞれに対応する第２の複数の評価値を取得する第１０のステップと、
    　前記第２の複数の評価値に基づいて、前記第２の複数の解の個数を超える個数の解を想定した場合の最大評価値を推定し、第２の最大評価値とする第１１のステップと、
    　前記第１の最適解の評価値を取得する第１２のステップと、
    前記第１の最適解の評価値が、前記第２の最大評価値を超えているか否かを判定する第１３のステップと、
    　を含む請求項１から４のいずれか１項に記載の最適解判定方法。
  6. 　前記第２の複数の解は、Ｐ，Ｑをそれぞれ自然数とすると、Ｐ×Ｑ個の解であり、
    　前記第１１のステップは、前記Ｐ×Ｑ個の解をＱ個のブロックに分け、前記ブロック毎のＰ個の解の評価値の区分最大値をＱ個取得し、前記Ｑ個の区分最大値を用いて、前記区分最大値が一般極値分布に従うものとして前記第２の最大評価値を推定する、
    　請求項５に記載の最適解判定方法。
  7. 　前記第１３のステップによる判定結果を出力する第１４のステップを更に含む、
    　請求項５又は６に記載の最適解判定方法。
  8. 　前記第１３のステップにより前記第１の最適解の評価値が前記第２の最大評価値を超えていないと判定されると、前記一定範囲を拡大し、前記拡大した一定範囲外の複数の解を第３の複数の解として一様抽出する第１５のステップと、
    　前記第３の複数の解のそれぞれに対応する第３の複数の評価値を取得する第１６のステップと、
    　前記第３の複数の評価値に基づいて、前記第３の複数の解の個数を超える個数の解を想定した場合の最大評価値を推定し、第３の最大評価値とする第１７のステップと、
    　前記第１の最適解の評価値が、前記第３の最大評価値を超えているか否かを判定する第１８のステップと、
    　を更に含む、
    　請求項５から７のいずれか１項に記載の最適解判定方法。
  9. 　前記第１３のステップにより前記第１の最適解の評価値が前記第２の最大評価値を超えていないと判定されると、前記解空間上の解であって、かつ前記第１の最適解から一定距離離れた第４の解候補を取得する第１９のステップと、
    　前記解空間上において、前記第１の最適解及び前記第４の解候補のそれぞれからの前記解空間上の距離が一定範囲外の第４の複数の解を一様抽出する第２０のステップと、
    　前記第４の複数の解のそれぞれに対応する第４の複数の評価値を取得する第２１のステップと、
    　前記第４の複数の評価値に基づいて、前記第４の複数の解の個数を超える個数の解を想定した場合の最大評価値を推定し、第４の最大評価値とする第２２のステップと、
    　前記第１の最適解の評価値が、前記第４の最大評価値を超えているか否かを判定する第２３のステップと、
    　を更に含む、
    　請求項５から７のいずれか１項に記載の最適解判定方法。
  10. 　前記解空間は、第１の制約条件における第１の解空間と第２の制約条件における第２の解空間とを含み、
    　前記第１の解空間に属する第１の解候補を取得する第２４のステップと、
    　前記第２の解空間上の複数の解を第５の複数の解として一様抽出する第２５のステップと、
    　前記第５の複数の解のそれぞれに対応する第５の複数の評価値を取得する第２６のステップと、
    　取得した前記第５の複数の評価値に基づいて、前記第５の複数の解の個数を超える個数の解を想定した場合の最大評価値を推定し、第５の最大評価値とする第２７のステップと、
    　前記第２の解空間における解候補であって、前記第２４のステップで取得した第１の解候補からの前記解空間上の距離が近い解を近傍解として取得する第２８のステップと、
    　前記近傍解の評価値を取得する第２９のステップと、
    　前記近傍解の評価値が前記第５の最大評価値の信頼区間内に入るか否かを判定する第３０のステップと、
    　を含む請求項１から９のいずれか１項に記載の最適解判定方法。
  11. 　前記第１のステップは、前記組合せ最適化問題における組合せ可能なパターンのうち、組合せの一部により残りの組合せを考慮せずとも不適となることが確定するかどうかを識別することで、識別すべきパターンを削減するステップ、及び前記組合せ可能なパターンのうち、組合せの一部だけに差分があるパターン群の共通部分を抽出し、残りの組合せを共有することで、識別すべきパターンを削減するステップのうちの少なくとも一方を用いて、前記組合せ可能パターンを縮約して列挙索引化するデータ構造を用い、前記解空間上の解の総数を算出し、前記算出した総数以下の乱数を発生させ、発生させた乱数により特定されるパターンに対応する解を抽出する、
    　請求項１から１０のいずれか１項に記載の最適解判定方法。
  12. 　組合せ最適化問題は、遺伝子制御ネットワークの組合せ最適化問題である請求項１から１１のいずれか１項に記載の最適解判定方法。
  13. 　請求項１から１２のいずれか１項に記載の最適解判定方法をコンピュータに実行させる最適解判定プログラム。
  14. 　請求項１３に記載の最適解判定プログラムを記録したコンピュータ読み取り可能な非一時的記録媒体。
  15. 　組合せ最適化問題における解の最適性を判定する最適解判定装置であって、
    　前記組合せ最適化問題の解空間上の複数の解を第１の複数の解として一様抽出する解抽出部と、
    　前記一様抽出した前記第１の複数の解のそれぞれに対応する評価値を取得する第１の評価値取得部と、
    　前記取得した前記第１の複数の評価値に基づいて、前記第１の複数の解の個数を超える個数の解を想定した場合の最大評価値を推定し、第１の最大評価値とする第１の最大評価値推定部と、
    　前記解空間に属する解のうち少なくとも１つの解を解候補として取得する解取得部と、
    　前記解候補に対応する評価値を前記第１の評価値取得部から取得し、前記取得した前記解候補に対応する評価値が前記第１の最大評価値の信頼区間内に入るか否かを判定し、前記解候補に対応する評価値が前記第１の最大評価値の信頼区間内に入ると判定された場合に、前記解候補を第１の最適解とする第１の判定部と、
    　を備えた最適解判定装置。
  16. 　前記解空間上の解であって、前記第１の最適解からの前記解空間上の距離が一定範囲外の複数の解を第２の複数の解として一様抽出する第２の解抽出部と、
    　前記第２の複数の解のそれぞれに対応する第２の複数の評価値を取得する第２の評価値取得部と、
    　前記第２の複数の評価値に基づいて、前記第２の複数の解の個数を超える個数の解を想定した場合の最大評価値を推定し、第２の最大評価値とする第２の最大評価値推定部と、
    　前記第１の最適解に対応する評価値を前記第２の評価値取得部から取得し、前記取得した前記第１の最適解に対応する評価値が、前記第２の最大評価値を超えているか否かを判定する第２の判定部と、
    　を更に備えた請求項１５に記載の最適解判定装置。
PCT/JP2018/006501 2017-03-15 2018-02-22 最適解判定方法、最適解判定プログラム及び最適解判定装置 WO2018168383A1 (ja)

Priority Applications (4)

Application Number Priority Date Filing Date Title
JP2019505816A JP6851460B2 (ja) 2017-03-15 2018-02-22 最適解判定方法、最適解判定プログラム、非一時的記録媒体及び最適解判定装置
CN201880009956.4A CN110249344A (zh) 2017-03-15 2018-02-22 最佳解判定方法、最佳解判定程序及最佳解判定装置
EP18768409.7A EP3598350A4 (en) 2017-03-15 2018-02-22 OPTIMAL SOLUTION EVALUATION METHOD, OPTIMAL SOLUTION EVALUATION PROGRAM, AND OPTIMAL SOLUTION EVALUATION DEVICE
US16/513,202 US11816580B2 (en) 2017-03-15 2019-07-16 Optimal solution determination method, optimal solution determination program, and optimal solution determination device

Applications Claiming Priority (2)

Application Number Priority Date Filing Date Title
JP2017050214 2017-03-15
JP2017-050214 2017-03-15

Related Child Applications (1)

Application Number Title Priority Date Filing Date
US16/513,202 Continuation US11816580B2 (en) 2017-03-15 2019-07-16 Optimal solution determination method, optimal solution determination program, and optimal solution determination device

Publications (1)

Publication Number Publication Date
WO2018168383A1 true WO2018168383A1 (ja) 2018-09-20

Family

ID=63522077

Family Applications (1)

Application Number Title Priority Date Filing Date
PCT/JP2018/006501 WO2018168383A1 (ja) 2017-03-15 2018-02-22 最適解判定方法、最適解判定プログラム及び最適解判定装置

Country Status (5)

Country Link
US (1) US11816580B2 (ja)
EP (1) EP3598350A4 (ja)
JP (1) JP6851460B2 (ja)
CN (1) CN110249344A (ja)
WO (1) WO2018168383A1 (ja)

Cited By (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2020201556A (ja) * 2019-06-06 2020-12-17 日本電信電話株式会社 近似ｚｄｄ構築方法、近似ｚｄｄ構築装置及びプログラム
EP3996011A4 (en) * 2019-07-04 2023-08-09 Daikin Industries, Ltd. COMBINATORY SOLUTION DETERMINATION SYSTEM

Families Citing this family (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20190155966A1 (en) * 2017-11-17 2019-05-23 Autodesk, Inc. Computer-implemented synthesis of a mechanical structure using a divergent search algorithm in conjunction with a convergent search algorithm
JP7283318B2 (ja) * 2019-09-12 2023-05-30 富士通株式会社 最適化装置、最適化プログラム、及び最適化方法

Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2000057124A (ja) * 1998-08-14 2000-02-25 Nec Corp 組合せ最適化方法および組合せ最適化システム
WO2004047040A1 (fr) 2002-11-13 2004-06-03 Sylvie Borne Systeme de securite pour des personnes susceptibles d'etre immergees
WO2004047020A1 (en) 2002-11-19 2004-06-03 Gni Usa Nonlinear modeling of gene networks from time series gene expression data

Family Cites Families (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US6865325B2 (en) * 2001-04-19 2005-03-08 International Business Machines Corporation Discrete pattern, apparatus, method, and program storage device for generating and implementing the discrete pattern
EP1436611A4 (en) * 2001-09-26 2007-04-11 Gni Kk BIOLOGICAL IDENTIFICATION USING GENREGULATION NETWORKS GENERATED FROM MULTIPLE INTERRUPTED EXTRESSION LIBRARIES

Patent Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2000057124A (ja) * 1998-08-14 2000-02-25 Nec Corp 組合せ最適化方法および組合せ最適化システム
WO2004047040A1 (fr) 2002-11-13 2004-06-03 Sylvie Borne Systeme de securite pour des personnes susceptibles d'etre immergees
WO2004047020A1 (en) 2002-11-19 2004-06-03 Gni Usa Nonlinear modeling of gene networks from time series gene expression data

Non-Patent Citations (6)

* Cited by examiner, † Cited by third party
Title
GIDDINGS, A.P.RARDIN, R.L.UZSOY, R.: "Statistical optimum estimation techniques for combinatorial optimization problems: a review and critique", JOURNAL OF HEURISTICS, vol. 20, no. 3, 2014, pages 329 - 358
GOLDEN, B.L.ALT, F.B.: "Interval estimation of a global optimum for large combinatorial problems", NAVAL RESEARCH LOGISTICS QUARTERLY, vol. 26, 1979, pages 69 - 77
IWASHITA, H.KAWAHARA, J.MINATO, S.: "TCS Technical Report, TCS-TR-A-12-60", 19 September 2012, HOKKAIDO UNIVERSITY, article "ZDD-Based Computation of the Number of Paths in a Graph"
NAKAGAWA RYOMA: "Performance Evaluation with Extreme Statistics on Approximate Solution of Traveling Salesman Problem", ABSTRACTS OF THE 2016 FALL NATIONAL CONFERENCE OF THE OPERATIONS RESEARCH SOCIETY OF JAPAN, THE 38TH CORPORATE CASE EXCHANGE CONFERENCE, September 2016 (2016-09-01), pages 174 - 175, XP009516462 *
SOMEYA HIROSHI: "Recent trends in development of stochastic optimization algorithms", THE PAPERS OF TECHNICAL MEETING ON SYSTEMS, IEE JAPAN, vol. ST-13, 27 June 2013 (2013-06-27), pages 51 - 55, XP009515654 *
TAGAWA, KIYOHARU: "Empirical distribution and solution by difference evolution of Optimisation problem including individual opportunity constraints", IPSJ SIG TECHNICAL REPORT, vol. 2017-MPS-112, 20 February 2017 (2017-02-20), pages 1 - 6, XP055609102, Retrieved from the Internet <URL:https://ipsj.ixsq.nii.ac.jp/ej/?action=repository_uri&item_id=177442&file_id=1&file_no=1> [retrieved on 20170314] *

Cited By (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2020201556A (ja) * 2019-06-06 2020-12-17 日本電信電話株式会社 近似ｚｄｄ構築方法、近似ｚｄｄ構築装置及びプログラム
JP7093972B2 (ja) 2019-06-06 2022-07-01 日本電信電話株式会社 近似ｚｄｄ構築方法、近似ｚｄｄ構築装置及びプログラム
EP3996011A4 (en) * 2019-07-04 2023-08-09 Daikin Industries, Ltd. COMBINATORY SOLUTION DETERMINATION SYSTEM

Also Published As

Publication number Publication date
JPWO2018168383A1 (ja) 2019-12-12
EP3598350A1 (en) 2020-01-22
US20190340513A1 (en) 2019-11-07
US11816580B2 (en) 2023-11-14
JP6851460B2 (ja) 2021-03-31
CN110249344A (zh) 2019-09-17
EP3598350A4 (en) 2020-04-22

Similar Documents

Publication Publication Date Title
Li et al. IsoLasso: a LASSO regression approach to RNA-Seq based transcriptome assembly
WO2018168383A1 (ja) 最適解判定方法、最適解判定プログラム及び最適解判定装置
Elmsallati et al. Global alignment of protein-protein interaction networks: A survey
Xu et al. Vfold: a web server for RNA structure and folding thermodynamics prediction
Feng et al. Inference of isoforms from short sequence reads
Maretty et al. Bayesian transcriptome assembly
JP6751376B2 (ja) 最適解探索方法、最適解探索プログラム及び最適解探索装置
González-Domínguez et al. GPU-accelerated exhaustive search for third-order epistatic interactions in case–control studies
Blazewicz et al. Graph algorithms for DNA sequencing–origins, current models and the future
JP7109339B2 (ja) ポリマー設計装置、プログラム、および方法
Nepomuceno et al. Pairwise gene GO-based measures for biclustering of high-dimensional expression data
Fonseca et al. Ranking beta sheet topologies with applications to protein structure prediction
Molinari et al. Transcriptome analysis using RNA-Seq fromexperiments with and without biological replicates: areview
Parikh et al. A data-driven architecture using natural language processing to improve phenotyping efficiency and accelerate genetic diagnoses of rare disorders
Liang et al. Endotyping in Heart Failure―Identifying Mechanistically Meaningful Subtypes of Disease―
Wu et al. TIGER: tiled iterative genome assembler
US20190205394A1 (en) Method for generating text string dictionary, method for searching text string dictionary, and system for processing text string dictionary
CN103310128A (zh) 考虑种子片段的长度的碱基序列处理系统及方法
Chen et al. Constructing consensus genetic maps in comparative analysis
JP2011024473A (ja) アプタマー分類装置、アプタマー分類方法、プログラムおよび記録媒体
McVicar et al. FPGA acceleration of short read alignment
Chang et al. An integer programming approach to DNA sequence assembly
Hua et al. PGS: a dynamic and automated population-based genome structure software
RU2799750C2 (ru) Биоинформационные системы, устройства и способы для выполнения вторичной и/или третичной обработки
Varma et al. Hardware acceleration of de novo genome assembly

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: 18768409

Country of ref document: EP

Kind code of ref document: A1

ENP Entry into the national phase

Ref document number: 2019505816

Country of ref document: JP

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: 2018768409

Country of ref document: EP

Effective date: 20191015