An efficient second-order cone programming approach for optimal selection in tree breeding

An important problem in tree breeding is optimal selection from candidate pedigree members to produce the highest performance in seed orchards, while conserving essential genetic diversity. The most beneficial members should contribute as much as possible, but such selection of orchard parents would reduce performance of the orchard progeny due to serious inbreeding. To avoid inbreeding, we should include a constraint on the numerator relationship matrix to keep a group coancestry under an appropriate threshold. Though an SDP (semidefinite programming) approach proposed by Pong-Wong and Woolliams gave an accurate optimal value, it required rather long computation time. In this paper, we propose an SOCP (second-order cone programming) approach to reduce this computation time. We demonstrate that the same solution is attained by the SOCP formulation, but requires much less time. Since a simple SOCP formulation is not much more efficient compared to the SDP approach, we exploit a sparsity structure of the numerator relationship matrix, and formulate the SOCP constraint using Henderson's algorithm. Numerical results show that the proposed SOCP approach reduced computation time in a case study from 39,200 seconds under the SDP approach to less than 2 seconds.


Introduction
The usage of mathematical optimization approaches for tree breeders have been increasing [1,15,20,23], since one of their purposes is to derive better performance from seed orchards. When tree breeders make a plan for new seed orchards, they determine the contributions of candidate pedigree members so that the resultant orchard maximizes response to the selection. From the viewpoint of mathematical optimization, the simplest form of the optimal selection problem is: Here, we assume that the number of the candidate members is m. The variable vector is x ∈ R m , and corresponds to the contributions of the candidates. The first constraint e T x = 1 indicates that the total contribution of candidate members is unity. We use the vector e ∈ R m to denote the vector of all ones, and the superscript T to denote the transpose of a vector or a matrix. In the second constraints, l ∈ R m and u ∈ R m are element-wise lower and upper bounds on the contributions, respectively. The performance measure appears in the objective function as its coefficient g = (g 1 , . . . , g m ) T , and the estimated breeding value (EBV) [13] is often employed for g. The values g 1 , . . . , g m are calculated separately before we solve the optimal selection problem, so we can consider g 1 , . . . , g m as constant values. The above problem is a simple linear programming problem, therefore, a greedy method is enough to solve it. Such solution includes the candidates corresponding to the highest EBV as much as possible, and it is most efficient in situations where all the pedigrees are independent. Even if the pedigrees are independent, diversity is still an issue. Lindgren et al. [11] discussed a linear deployment in which the contributions of the candidate members are proportional to their EBVs. In practical situations, however, we cannot overlook the effects due to the relatedness that accumulates over cycles of breeding the pedigrees.
To reflect the effect of the relatedness, Cockerham [4] extended the definitions of coancestry coefficients in order to include coancestry of a group. The group coancestry on the contributions x is calculated with the formula x T Ax 2 , where A ∈ R m×m is the numerator relationship matrix of Wright [29] (we will review a formula for the numerator relationship matrix in Section 2). Introducing a constraint to keep group coancestry under an appropriate level θ ∈ R, Meuwissen [15] proposed a formulation of optimal contributions: max : g T x subject to : Meuwissen developed an iterative method based on Lagrangian multipliers to solve this optimization problem, and his method has been used in breedings, for example, [6,9,28]. It is a characteristic of this method that some variables x i may be fixed to its lower or upper bounds (l i or u i ) during the iterations, and Pong-Wong and Woolliams [20] demonstrated that the method did not always obtain the optimal solution.
Pong-Wong and Woolliams utilized the structure of the numerator relationship matrix A to formulate the problem (1) into a semidefinite programming (SDP) problem. SDP is a convex optimization problem that maximizes a linear objective function over the constraints described as linear matrix inequalities. The research in 1990s, for example [7,10], extended interior-point methods from linear programming problems to SDP problems. Based on primal-dual interiorpoint methods, software packages (like SDPA [30], SDPARA [31], SDPA-C [32], SDPT3 [26], and SeDuMi [25]) have been developed for solving SDPs. Using the SDP formulation, Pong-Wong and Woolliams obtained the exact optimal value of (1). The number of candidate members discussed in [20], however, was limited to only small sizes, m ≤ 10. Ahlinder et al. [1] implemented the SDP approach into a software package called OPSEL [17] 4 with the help of the latest version of SDPA (a high-performance SDP solver) [30]. They solved large problems (m ≥ 10, 000) that were generated from real Scots pine pedigrees and performance data, and they also focused on flexibility and re-optimization of the SDP approach.
A main obstacle in the numerical tests of [1] was that the SDP approach took rather long computation time. They reported for their case study that it required five-hours of computation time for a problem of the size m = 12, 000. Even though the SDP guarantee the optimal solution, the computation time is rather long for operational application and requiring significant truncation of the candidate list prior to optimizing the selection.
In this paper, we propose a second-order cone programming (SOCP) approach. SOCP is a convex optimization that maximizes a linear objective function over second-order cone constraints. SOCP can be considered as a special case of SDP, and can be efficiently solved with interior-point methods in a similar way to SDP [24,27]. Lobo et al. [12] discussed wide-range applications of SOCP, for example, filter design and truss design, and Sasakawa and Tsuchiya applied SOCP to magnetic shield design [22]. Alizadeh and Goldfarb [2] surveyed theoretical and algorithmic aspects of SOCP, and the software packages SDPT3 [26] and SeDuMi [25] can solve not only SDP but also SOCP using the primal-dual interior-point methods. In addition, ECOS [5] was also implemented recently to solve SOCP problems.
We first discuss that the proposed SOCP formulation also attains the optimal solution of (1). Since SOCP is a special case of SDP, we could expect that a simple SOCP formulation would be enough to reduce the computation time. However, preliminary numerical tests showed that the simple formulation did not perform well. We therefore utilize a sparsity embedded in the numerator relationship matrix and establish a more efficient SOCP formulation. Furthermore, we integrate Henderson's algorithm [8] into this formulation. Numerical tests with the data including Scots pine showed that the SOCP formulation with Henderson's algorithm reduced a great amount of computation time. For the case of m = 10, 100, we attained a speedup of 20,000-times compared to the SDP approach.
The rest of this paper is organized as follows. Section 2 describes the SDP approach of Pong-Wong and Woolliams and discusses a simple SOCP formulation. In Section 3, we propose SOCP formulations and derive an efficient method to solve problem (1) . Section 4 shows the numerical results to verify the computation time reduction for problems of various sizes. Finally, Section 5 gives conclusions and discusses future directions.

SDP formulation and simple SOCP formulation
Since a principal characteristic of our problem (1) is determined by the numerator relationship matrix A, we first review a formula to evaluate its elements. We then describe the SDP approach of Pong-Wong and Woolliams [20], and compare the performance of the SDP approach and a simple SOCP formulation.
To evaluate the elements of the numerator relationship matrix A, we separate the set of pedigree candidate members P := {1, 2, . . . , m} into the three disjoint groups: where    P 0 = {i ∈ P : both parents p(i) and q(i) are unknown} P 1 = {i ∈ P : one parent p(i) is known and the other parent q(i) is unknown} P 2 = {i ∈ P : both parents p(i) and q(i) are known}. Figure 1 gives an example of pedigree with m = 9 members and illustrates its genealogical chart. In this example, P 0 = {1, 2}, P 1 = {5}, P 2 = {3, 4, 6, 7, 8, 9}, and the parents of the 8th member are p(8) = 7 and q(8) = 6. We use a convention p(i) = 0 or q(i) = 0 if the parent p(i) or q(i) is unknown, respectively, and we can assume i > p(i) ≥ q(i) for all i ∈ P without loss of generality.
where we use a convention A pq = 0 if p = 0 or q = 0. When we apply this calculation to the example of Figure 1, we obtain the corresponding matrix A as follow: Pong-Wong and Woolliams [20] formulated the problem (1) into an SDP problem. A standard form of SDP is given by max : We use S n to denote the space of n × n symmetric matrices, and S n + ⊂ S n to denote the space of positive semidefinite matrices of dimension n. The variables are z 1 , . . . , z m ∈ R, and the input data are c 1 , . . . , c m ∈ R and F 0 , The key step of Pong-Wong and Woolliams [20] was the usage of the Schur complement of a matrix block. They noticed that the numerator relationship matrix A is always positive definite, and they utilized this property to convert the constraint on the group coancestry to a positive semidefinite condition on the symmetric matrix: With this positive semidefinite constraint, they converted (1) into the standard SDP form. The input vector c and the input matrices F 0 , F 1 , . . . , F m in (3) are given as follow [1,20]: We use e i to denote the vector of all zeros except 1 at the ith element, and Diag(u) to denote the diagonal matrix whose diagonal elements are u. Note that the dimension of F 0 , F 1 , . . . , F m is 3m + 3. The software package OPSEL [17] automatically formulates the optimal selection into an SDP problem of this form. Table 1 shows the computed optimal values and the computation time of Meuwissen's implementation (GENCONT) [16] and the SDP formulation. We executed the numerical tests using Matlab R2015a on Windows 8.1 PC with Xeon CPU E3-1231 (3.40 GHz, 4 cores) and 8 GB memory space. We used Windows, since GENCONT can run only on Windows. To solve the SDP (3), we employed SDPA [30].
We observe from Table 1 that the SDP formulation attained a better optimal value than GEN-CONT. Actually, as shown in [20], the optimal value of the SDP formulation was the exact optimal value, while the Lagrangian multiplier method implemented in GENCONT could not guarantee the optimality. For the large problem (m = 10, 100), GENCONT gave up the computation, but the SDP formulation again obtained the optimal solution. On the other hand, the disadvantage of the SDP formulation is its computation time. Even using the parallel computing implemented in SDPA, the SDP formulation was slower than GENCONT for m = 2, 045. Furthermore, the computation time for the large problem exceeded 10 hours even with four cores. When we used only one core, the SDP formulation would require longer computation time than 24 hours.
To reduce the heavy computation time of the SDP formulation, we review the group coancestry constraint x T Ax 2 ≤ θ. With the positive definiteness of A, we focus on the property that this constraint can also be described as a second-order condition. The vector v ∈ R nq is said to satisfy the second-order condition if v ∈ K nq . The symbol K nq denotes the second-order cone of dimension n q : The second-order cone is a special case of positive semidefinite constraint. In fact, using the Schur complement, we can verify that Furthermore, since the numerator relationship matrix A is positive definite, we can apply the Cholesky factorization to A to obtain the upper triangular matrix U such that A = U T U . This factorization derives a second-order cone condition that corresponds to the group coancestry constraint: Since the most difficult constraint in (1) can be expressed using a second-order cone, it is natural to consider its SOCP formulation. A standard form of SOCP problem is described by We use R n l + × K nq to denote the Cartesian product of the non-negative orthant of dimension n l , R n l + := {v ∈ R n l : v i ≥ 0 for i = 1, . . . , n l }, and the second-order cone of dimension n q . In this problem, the variable vector is z ∈ R m , and the input data are c ∈ R m , f 0 ∈ R n l +nq + and F ∈ R (n l +nq)×m . A more general SOCP form than (6) can be defined so that it can handle the Cartesian product of multiple second-order cones, but one cone is enough for the discussion in this paper.
If we formulate the optimal selection problem (1) in a simple style, we obtain an SOCP problem: We use I to denote the identity matrix of dimension m. We call this formulation a simple SOCP formulation. We emphasize that this simple SOCP formulation also gives the exact optimal value of the optimal selection problem (1) due to the equivalence from (4) and (5); Since SOCP is a special case of SDP, we expected that we could solve the SOCP formulation (7) faster than the SDP formulation (3). In Table 2, we compare the computation times of the SDP and SOCP formulations. We used ECOS [5] as the SOCP solver for (7). For the small problem (m = 2, 045), the SOCP formulation successfully reduced the computation time from 70.21 seconds to 0.28 seconds, so its speed-up was 250-times. However, for the large problem (m = 10, 100), the speed-up was limited to 7-times. When we consider the computational complexity, the simple SOCP formulation would be slower than the SDP formulation for further large problems. On contrary to the fact that SOCP is a special case of SDP, the simple SOCP formulation did not seem promising.

Efficient formulation based on SOCP
To obtain an efficient formulation based on SOCP, we investigated key properties of the simple SOCP formulation (7). In particular, we focused on the structure of the numerator relationship matrix A and its inverse A −1 , since the SDP formulation uses only A −1 . Figure 2 illustrates the positions of the non-zero elements of A and A −1 for the problem of size m = 10, 100. The dimensions of A and A −1 corresponds to the size m. We can see from Figure 2 that A −1 is much sparser than A. The numbers of non-zero elements in A and A −1 are 56,754,980 and 56,092, respectively, hence their density against the fully-dense matrix (m 2 non-zero elements) are 55.6% and 0.0549%, respectively. When we apply the Cholesky factorization to A, the upper-triangular matrix U inherits the dense property. The number of non-zero elements in U is 23,171,296, and its density against the fully-dense upper triangular matrix is 45.4%. We should utilize the property that A −1 is remarkably sparse compared to A and U . From this observation, the direction we should pursue is to use A −1 instead of A itself. A key step of our approach is to introduce the new variable y = Ax. Then, the replacement x by A −1 y in the optimal selection (1) leads to an equivalent optimization problem: The Cholesky factor of A −1 is the transposed matrix of U −1 , denoted by U −T : namely, A −1 = (U −T ) T U −T . In practical implementation, since A −1 is remarkably sparse, we apply an appropriate row/column permutation like AMD (approximate minimum degree permutation) [3] to A −1 in order to reduce the number of fill-in (the non-zero elements that newly appear in U −T during the process of the Cholesky factorization). Using U −T , we convert this problem into an SOCP We call this formulation a sparse SOCP formulation of the optimal selection problem (1). We use the matrix A to formulate this new SOCP formulation, but we do not have to use A to solve the resultant SOCP problem with the interior-point methods. Furthermore, when we obtain the optimal solution y * of (8), we can also obtain the optimal solution x * of the original problem (1) via the relation x * = A −1 y * without using the dense matrix A.
In Table 3, we compare the performance of the SDP formulation (3), the simple SOCP formulation (7), and the sparse SOCP formulation (8). In this table, the first column 'nnz' is the total number of non-zero elements of F 0 , F 1 , . . . , F m for SDP and that of f 0 and F for SOCP. The second column is the computation time to convert the pedigree like Figure 1 and EBVs (estimated breeding values) into the SDP or SOCP formulations, and the third column is the computation time of the solvers. We applied SDPA with four cores to the SDP formulation and ECOS to the SOCP formulations. The fourth (last) column is the total computation time. Table 3 shows that for the large problem, the computation time of the sparse SOCP formulation was much shorter than that of the simple SOCP formulation. The sparse SOCP formulation seemed to have more complex structure than the simple SOCP formulation since it repeatedly contained A −1 in the matrix F , but the total number of non-zeros was reduced from 23,231,801 in the simple SOCP formulation to 159,570 in the sparse SOCP formulation. This led the computation time reduction for the SOCP solver.
To reduce the total time further, we now address the computation time to build the SOCP formulation. In the case of m = 10, 100, the conversion time occupied 94 % of the total time. In particular, the construction of the dense matrix A and its inversion are the principal bottlenecks.
We investigate further the properties of A −1 , and employ a compact algorithm to construct A −1 proposed by Henderson [8]. In the compact algorithm, we use the vector of inbreeding   [21] devised an efficient method to compute the inbreeding coefficients h without constructing the matrix A itself, and Masuda et al. utilized this method to implement their YAMS package [14]. The compact algorithm in [8] with the enhancement [21] is summarized in Table 4. We use an indicator function δ(p) such that δ(0) = 1, and δ(p) = 0 for p = 0. We also use the notation A −1 ij to denote the (i, j) element of A −1 . The compact formula for Table 4 can be described as When we apply this formula to the example in Figure 1, we obtain the inverse of the numerator relationship matrix as follow: We can see in this example that A −1 has more zero-elements than A of (2). It is known that b i > 0 for any i = 1, 2, . . . , m from properties of the inbreeding coefficients. Therefore, the constraint on the group coancestry can be transformed into another second-order cone constraint: Here, B is the matrix whose ith row vector is defined by We now replace the matrix U −T in the sparse SOCP formulation (8) by B, hence, we derive another SOCP formulation: We call this formulation a compact SOCP formulation.
The combination of the compact algorithm in Table 4 and the second-order cone constraint (9) gives us a formulation that does not rely on any dense matrices. Table 5 adds the results of the compact SOCP formulation to Table 3. From Table 5, we observe that the solver time for the compact SOCP formulation was slightly shorter than the sparse SOCP formulation. A further computation time reduction was obtained in the computation time to build the SOCP formulations from the pedigree. In the case m = 10, 100, the conversion was reduced from 24.14 seconds to 0.37 seconds. Furthermore, we saved a lot of memory space. For the case m = 10, 100, the matrix A required 778 MB of memory, while in contrast the memory required for B is only 2.33 MB.

Numerical tests
We conducted a numerical evaluation of the SDP and SOCP formulations on several datasets. The datasets are for problems of sizes 2045, 5050, 15100, 15222, 50100, 100100 and 300100. The data with the sizes 2,045 and 15,222 were from Scots pine orchards and loblolly pine orchards, respectively, and these data are available at the Dryad Digital Repository http://dx.doi.org/ 10.5061/dryad.9pn5m. The other data were generation by simulation of five cycles of breeding in a closed population using the approach of [18,19].
Although we used Matlab for the comparison between the formulations, we also implemented the compact SOCP formulation with C++. There were two reasons to implement it outside a Matlab environment. The first is that we can directly know the structure of non-zero elements that appear in the matrix B. Therefore, a specified data structure can accelerate the computation time to arrange the input data for building F . Secondly, we expect a software package that does not depend on commercial software would extend the opportunity for the field application in tree breeding. Due to the latter reason, we also excluded commercial SOCP solvers from the numerical tests. One of the advantages of the compact SOCP formulation is that we do not need numerical routines for the Cholesky factorization. If we implement the simple or sparse SOCP formulation with C++, we have to embed certain numerical routines for the Cholesky factorization. In addition, the sparse Cholesky factorization requires a preprocessing by AMD [3] to derive its best performance. In contrast, the compact SOCP formulation obtains the matrix directly F as discussed in Section 3.
The computation environment for the small problems (m ≤ 10, 100) was Matlab R2015a on a Windows 8.1 PC with Xeon E3-1231 (3.40 GHz) and 8 GB memory space. For the large problems (m ≥ 15, 100), we used Matlab R2014b on a Debian Linux server with Opteron 4386 (3.10 GHz) and 128 GB memory space, since 8 GB memory space was not enough for the SDP formulation. Table 6 shows the numerical results of the SDP and SOCP formulations. We observe from this table that the SDP formulation demanded rather long computation time, particularly for the larger problems. For the problem m = 15, 222, the compact SOCP formulation with C++ reduced the 22, 566 seconds of the SDP formulation to only 2.21 seconds.
Another significant advantage of the compact formulation is memory consumption. The SDP formulation consumed 31 GB memory space to solve the problem m = 15, 222, and it failed to handle m ≥ 50, 100, despite having 128 GB of memory available. If we use a rough estimation, the memory required for the largest problem m = 300, 100 would be 12,000 GB. The simple SOCP formulation also suffered from a heavy memory requirement to store the dense matrix U . The sparse SOCP formulation did not require the dense matrix A in the resultant SOCP problem, but it utilized A and computed A −1 , hence it required the long conversion time in the same way as the SDP formulation. In contrast, the compact SOCP formulation consumed less than 766 MB memory space to solve even the largest problem m = 300, 100. This memory reduction was mainly a result of employing Henderson's algorithm.
From the numerical results, we also observe that the compact SOCP formulation with C++ is faster than that with Matlab. The discrepancy between Matlab (99.65 seconds) and C++ (82.41 seconds) in the largest problem was due to a specified data structure written in C++. In particular, the structure was effective when we arranged the pedigree before building F . * OOM = "out of memory" * > 24 hours = "the computation failed to complete with in in 24 hours."

Conclusions and future directions
We examined the SOCP formulations for the optimal selection problem arising from tree breeding. We employed the transformation x = A −1 y based on the sparsity of A −1 and the efficient method to build A −1 by the Henderson's algorithm with the Quaas enhancement. The compact SOCP formulation thus did not involve any dense matrix or the Cholesky factorization. The numerical results demonstrated that the compact SOCP formulation obtained the optimal solution significantly faster than the existing SDP formulation. The SOCP formulations proposed in this paper may look rather simple for researchers in the field of mathematical optimization. However, the computation time reduction in the optimal selection problem will help improve operational application in tree breeding. We expect that this paper will be one of bridges to introduce efficient approaches cultivated in mathematical optimization to tree breeders.
In this paper, we discussed an unequal deployment of parental genotypes to seed orchards, where the contributions of selected members are not required to be equal. To deal with equal deployment, as might be appropriate for selection of a fixed-size breeding population, we will need a method to solve a mixed integer SOCP problem; that is an SOCP problem in which some variables are constrained to be integers. The structure of the SOCP formulation developed in this paper will be a basis for an efficient method to consider the mixed integer SOCP problem in the optimal selection problems.