next up previous
Next: The primal-dual path-following algorithm Up: Infeasible-interior-point algorithms Previous: The search direction

Computation of specific search directions

In this package, the user has four choices of symmetrizations resulting in four different search directions, namely,

(1) the AHO direction, corresponding to P=I;
(2) the HKM direction, corresponding to tex2html_wrap_inline2205 ;
(3) the NT direction, corresponding to tex2html_wrap_inline2207 , where N is the unique matrix such that tex2html_wrap_inline2211 , and D is a diagonal matrix; and
(4) the GT direction, corresponding to tex2html_wrap_inline2215 , where the matrices D and tex2html_wrap_inline2219 are defined as follows. Suppose tex2html_wrap_inline2221 and tex2html_wrap_inline2223 are the Cholesky factorizations of X and Z respectively. Let the SVD of tex2html_wrap_inline2229 be tex2html_wrap_inline2231 . Let tex2html_wrap_inline2233 and tex2html_wrap_inline2235 be positive diagonal matrices such that the equalities tex2html_wrap_inline2237 and tex2html_wrap_inline2239 hold, with all the rows of tex2html_wrap_inline2219 and tex2html_wrap_inline2243 having unit norms. Then tex2html_wrap_inline2245 .

To describe our implementation 3, a discussion on the efficient computation of the Schur complement matrix M is necessary, since this is the most expensive step in each iteration of our algorithms where usually at least 80% of the total CPU time is spent. From equation (10), it is easily shown that the (i,j) element of M is given by


Thus for a fixed j, computing first the matrix tex2html_wrap_inline2257 and then taking inner product with each tex2html_wrap_inline2259 , tex2html_wrap_inline2261 , give the jth column of M.

However, the computation of M for the four search directions mentioned above can also be arranged in a different way. The operator tex2html_wrap_inline2269 corresponding to these four directions can be decomposed generically as


where tex2html_wrap_inline2271 denotes the Hadamard (elementwise) product and the matrices R, T, tex2html_wrap_inline2277 , and tex2html_wrap_inline2279 depend only on X and Z. Thus the (i,j) element of M in (14) can be written equivalently as


Therefore the Schur complement matrix M can also be formed by first computing and storing tex2html_wrap_inline2291 for each tex2html_wrap_inline2293 , and then taking inner products as in (15).

Computing M via different formulas, (14) or (15), will result in different computational complexities. Roughly speaking, if most of the matrices tex2html_wrap_inline2161 are dense, then it is cheaper to use (15). However, if most of the matrices tex2html_wrap_inline2161 are sparse, then using (14) will be cheaper because the sparsity of the matrices tex2html_wrap_inline2161 can be exploited in (14) when taking inner products. For the sake of completeness, in Table 1, we give an upper bound on the complexity of computing M for the above mentioned search directions when computed via (14) and (15).

Table 1: Upper bounds on the complexities of computing M (for real SDP data) for various search directions. We count one addition and one multiplication each as one flop. Note that, except for the HKM direction, all the others require an eigenvalue decomposition of a symmetric matrix in the computation of M.

The reader is referred to [2], [8], and [9] for more computational details, such as the actual formation of M and tex2html_wrap_inline2331 , involved in computing the above search directions. The derivation of the upper bounds on the computational complexities of M computed via (14) and (15) is given in [9]. The issue of exploiting the sparsity of the matrices tex2html_wrap_inline2161 is discussed in full detail in [5] for the HKM and NT directions; whereas an analogous discussion for the AHO and GT directions can be found in [9].

In our implementation, we decide on the formula to use for computing M based on the CPU time taken during the third and fourth iteration to compute M via (15) and (14), respectively. We do not base our decision on the first two iterations for two reasons. Firstly, if the initial iterates tex2html_wrap_inline2341 and tex2html_wrap_inline2343 are diagonal matrices, then the CPU time taken to compute M during these two iterations would not be accurate estimate of the time required for subsequent iterations. Secondly, there are overheads incurred when variables are first loaded into MATLAB workspace.

We should mention that the issue of exploiting the sparsity of the matrices tex2html_wrap_inline2161 is not fully resolved in our package. What we have used in our package is only a rough guide. Further improvements on this topic are left for future work.

next up previous
Next: The primal-dual path-following algorithm Up: Infeasible-interior-point algorithms Previous: The search direction