Reflection methods for inverse problems with application to protein conformation determination
CChapter 1
Reflection methods for inverse problems withapplications to protein conformationdetermination
Jonathan M. Borwein and Matthew K. Tam
Abstract
The Douglas–Rachford reflection method is a general purpose algorithmuseful for solving the feasibility problem of finding a point in the intersection offinitely many sets. In this chapter we demonstrate that applied to a specific prob-lem, the method can benefit from heuristics specific to said problem which exploitits special structure. In particular, we focus on the problem of protein conforma-tion determination formulated within the framework of matrix completion, as wasconsidered in a recent paper of the present authors.
Key words: reflection methods; inverse problems; protein conformation
This chapter builds on a series of seven lectures titled
Techniques of VariationalAnalysis given by the first author at the CIMPA school
Generalized Nash Equilib-rium Problems, Bilevel Programming and MPEC held November 25 to December6, 2013, University of Delhi, New Delhi, India. In this written presentation we focuson reflection methods for protein conformation determination , as was discussed inthe seventh and final lecture of the series. The complete lectures — one through sixtaken from [13] — can be found online at:
Before turning our attention to reflection methods, we briefly outline the contentof the first six lectures.
J.M. BorweinCARMA Centre, University of Newcastle, Callaghan, NSW 2308, Australia. e-mail: [email protected]
M.K. TamCARMA Centre, University of Newcastle, Callaghan, NSW 2308, Australia. e-mail: [email protected] a r X i v : . [ m a t h . O C ] A ug J.M. Borwein, M.K. Tam • Lectures 1 & 2 provided an introduction to variational analysis and variationalprinciples [13, § § • Lectures 3 & 4 introduced nonsmooth analysis: normal cones and subdifferen-tials of lower semi-continuous functions, Fr´echet and limiting calculus [13, § § § § • Lecture 5 turned to multifunction analysis: sequences of sets, continuity of maps,minimality and maximal monotonicity, and distance functions [13, § § • Lecture 6 focussed on convex feasibility problems and the method of alternatingprojections [13, § Given a (finite) family of sets, the corresponding feasibility problem is to find a pointcontained in their intersection.
Douglas–Rachford reflection methods form a classof general-purpose iterative algorithms which are useful for solving such problems.At each iteration, these methods perform (metric) reflections and (metric/nearestpoint) projections with respect to the individual constraint sets in a prescribed fash-ion. Such methods are most useful when applied to feasibility problems whose con-straint sets have more easily computable reflections and projections than does theintersection.When the underlying constraint sets are all convex, Douglas–Rachford methodsare relatively well understood [6, 12, 11, 7] — their behaviour can be analysed usingnonexpansivity properties of convex projections and reflections. In the absence ofconvexity, recent result have assumed the constraint sets to possess other structuraland regularity properties [10, 1, 20]. However, at present, this developing theoreticalfoundation is not sufficiently rich to explain many of the successful applications inwhich one or more of the constraint sets lacks convexity [3, 2, 17, 18]. In thesecases, the method can be viewed as a heuristic inspired by its behaviour within fullyconvex settings.More generally, with any algorithm there is typically a trade-off between thescope of their applicability and tailoring of performance to particular instances.Douglas–Rachford reflection methods are no different. Owing to these methods’broad applicability, potential for further problem specific refinements when appliedto special classes of feasibility problems is possible.In this chapter, we investigate and develop one such refinement with a focus onapplication of the Douglas–Rachford method to protein conformation determina-tion . This application was previously considered as part of [3]. We now proposeproblem specific heuristics, and also study the effect of increasing problem size. Wefinish by demonstrating a complementary application of the approach arising in thecontext of ionic liquid chemistry .The remainder of this chapter is organized as follows. In Sections 1.3, 1.4, 1.5 &1.6 we introduce the necessary mathematical preliminaries along with the Douglas–
Reflection methods for inverse problems 3
Rachford reflection method, before formulating the protein conformation determi-nation problem. Substantial numerical and graphical results are given in Section 1.7,and concluding remarks in Section 1.8.
Let E denote a Euclidean space, that is, a finite dimensional Hilbert space. We willmainly be concerned with the space R m × m (i.e., real m × m matrices) equipped withthe inner-product given by (cid:104) A , B (cid:105) : = tr ( A T B ) . Here the symbol tr ( X ) (resp. X T ) denotes the trace (resp. transpose) of the matrix X . The induced norm is the Frobenius norm and can be expressed as (cid:107) A (cid:107) F : = (cid:113) tr ( A T A ) = (cid:115) m ∑ i = m ∑ j = a i j . The subspace of real symmetric m × m matrices is denoted S m , and the cone ofpositive semi-definite m × m matrices by S m + .Given sets C , C , . . . , C N ⊆ E , the feasibility problem isfind x ∈ N (cid:92) i = C i . (1.1)When the intersection in (1.1) is empty, one often seeks a “good” surrogate for apoint in the intersection. When N =
2, a useful surrogate is a pair of points, onefrom each set, which minimize the distance between the sets – a best approximationpair [6]. A partial (real) matrix is an m × m array for which entries only in certain locationsare known. Given a partial matrix A = ( a i j ) ∈ R m × m , a matrix B = ( b i j ) ∈ R m × m isa completion of A if b i j = a i j whenever a i j is known. The problem of (real) matrixcompletion is the following: Given a partial matrix find a completion belonging toa specified family of matrices .Matrix completion can be naturally formulated as a feasibility problem. Let A bethe partial matrix to be completed. Choose C , C , . . . , C N such that their intersectionis equal to the intersection of completions of A with the specified matrix family.Then (1.1) is precisely the problem of matrix completion for A . The simplest such J.M. Borwein, M.K. Tam case is when C is the set of all completions of A and the intersection of C , . . . , C N equals the desired matrix class. Remark 1.
More generally, one may profitably consider matrix completion for rect-angular matrices [3], for example with doubly stochastic matrices. However, sincethe partial matrices in the discussed protein application are always square, for thepurposes of this discussion, we only concern ourselves with the square case.
The projection onto C ⊆ E is the set-valued mapping P C : E ⇒ C which maps anypoint x ∈ E to its sets of nearest points in C . More precisely, P C ( x ) = (cid:26) c ∈ C : (cid:107) x − c (cid:107) ≤ inf y ∈ C (cid:107) x − y (cid:107) (cid:27) . The reflection with respect to C is the set-valued mapping R C : E ⇒ E given by R C = P C − I , where I denotes the identity mapping.When C is non-empty, closed, and convex, its corresponding projection operator(and hence its reflection) is single-valued (see, for example, [15, Ch. 1.2]). x p r x p r x p p r r Fig. 1.1 (Left) The (single-valued) projection, p i , and reflection, r i , of the point x i onto a convexset, for i = ,
2. (Right) The (set-valued) projection, { p , p } , and reflection, { r , r } , of the point x onto a non-convex set. Note the non-expansivity of the reflection in the convex case. Given A , B ⊆ E and x ∈ E , the Douglas–Rachford reflection method is the fixedpoint iteration given by x n + ∈ T A , B x n where T A , B = I + R B R A . (1.2)We refer to the sequence ( x n ) ∞ n = as a Douglas–Rachford sequence , and to the map-ping T A , B as the Douglas–Rachford operator .We now recall the behavior of the Douglas–Rachford method in the classicalconvex setting. In this case, T A , B is single-valued as a consequence of the single- Reflection methods for inverse problems 5 valuedness of each of P A , P B , R A and R B . We denote the set of fixed points of asingle-valued mapping T by Fix T = { x ∈ E : T x = x } , and the normal cone of aconvex set C at the point x by N C ( x ) = (cid:40) { y ∈ E : (cid:104) C − x , u (cid:105) ≤ } if x ∈ C , /0 otherwise.For convenience, we also introduce the two sets E = (cid:26) x ∈ A : inf a ∈ A (cid:107) a − x (cid:107) ≤ inf a ∈ A , b ∈ B (cid:107) a − b (cid:107) (cid:27) , F = (cid:26) x ∈ B : inf b ∈ B (cid:107) x − b (cid:107) ≤ inf a ∈ A , b ∈ B (cid:107) a − b (cid:107) (cid:27) , and the vector v = P B − A ( ) . Here the overline denotes the closure of the set. Theorem 1 (Convex Douglas–Rachford in finite dimensions [6]).
Suppose A , B ⊆ E are closed and convex. For any x ∈ E define x n + = T A , B x n . Then there is some v ∈ E such that:(i) x n + − x n = P B R A x n − P A x n → v and P B P A x n − P A x n → v.(ii) If A ∩ B (cid:54) = /0 then ( x n ) ∞ n = converges to a point in Fix ( T A , B ) = ( A ∩ B ) + N A − B ( ) ; otherwise, (cid:107) x n (cid:107) → + ∞ .(iii) Exactly one of the following two alternatives holds.(a) E = /0 , (cid:107) P A x n (cid:107) → + ∞ , and (cid:107) P B P A x n (cid:107) → + ∞ .(b) E (cid:54) = /0 , the sequences ( P A x n ) ∞ n = and ( P B P A x n ) ∞ n = are bounded, and theircluster points belong to E and F, respectively; in fact, the cluster points of (( P A x n , P B R A x n )) ∞ n = and (( P A x n , P B P A x n )) ∞ n = are a best approximation pairs relative to ( A , B ) . Theorem 1 provides the template for application of the Douglas–Rachford methodas a heuristic for non-convex feasibility problems. Furthermore, this theorem alsoshows that for the Douglas–Rachford method the sequence of primary interest is notthe fixed point iterates ( x n ) ∞ n = themselves, but their shadows ( P A x n ) ∞ n = . Remark 2 (Douglas–Rachford splitting).
The Douglas–Rachford reflection methodcan be viewed as a special case of the
Douglas–Rachford splitting algorithm forfinding a zero of the sum of two maximally monotone operators. This more generalsplitting method iterates by using resolvents of the given maximally monotone op-erators rather than projection operators of sets. The reflection method is obtained inthe special case in which the maximal monotone operators are normal cones to thefeasibility problem sets. For details, we refer the reader to [5].
J.M. Borwein, M.K. Tam x n R A x n R B R A x n x n + = T A , B x n AB Fig. 1.2
One iteration of the Douglas–Rachford method for the sets A = { x ∈ E : (cid:107) x (cid:107) ≤ } and B = { x ∈ E : (cid:104) a , x (cid:105) = b } . Within an implementation of the Douglas–Rachford method, computation of theprojection operators are the component typically requiring the most resources. Itis therefore beneficial to store two additional sequences in memory; the shadowsequence ( P A x n ) ∞ n = , and the sequence ( P B R A x n ) ∞ n = . This is because iteration (1.2)is expressible as x n + ∈ x n + P B R A x n − P A x n = x n + P B ( P A x n − x n ) − P A x n . (1.3)An implementation utilizing this approach is given in Algorithm 1.3. The stoppingcriterion uses a relative error and is discussed in Section 1.7. Fig. 1.3
Implementation of the basic Douglas–Rachford algorithm.
Input : x ∈ E and ε > n = p ∈ P A ( x ) ; while n = or (cid:107) r n − p n (cid:107) > ε (cid:107) p n (cid:107) do r n ∈ P B ( p n − x n ) ; x n + = x n + r n − p n ; p n + ∈ P A ( x n + ) ; n = n + endOutput : p n ∈ E Reflection methods for inverse problems 7
Proteins are large biomolecules which are comprised of multiple amino acid residues, each of which typically consists of between 10 and 25 atoms. Proteins participate isvirtually every cellar process, and knowledge of their structural conformation givesinsight into the mechanisms by which they perform.One of many techniques that can be used to determine conformation is nuclearmagnetic resonance (NMR) . Currently NMR is only able to non-destructively re-solve relatively short distances ( i.e., those less than ∼ low-rank Euclidean distance matrix completion problem . We next introduce the neces-sary definitions.We say that a matrix D = ( D i j ) ∈ R m × m is a Euclidean distance matrix (EDM) ifthere exists points z , z , . . . , z n ∈ R m such that D i j = (cid:107) z i − z j (cid:107) for i , j = , , . . . , m . (1.4)Clearly any EDM is symmetric, non-negative, and hollow ( i.e., contains only zerosalong its main diagonal). When (1.4) holds for a set of points in R q , we say D is embeddable in R q . If D is embeddable in R q but not in R q − , then we say that D is irreducibly embeddable in R q .We now recall a useful characterization of EDMs, due to Hayden and Wells [19].In what follows, the matrix Q ∈ R m × m is the Householder matrix given by Q = I − vv T v T v , where v = (cid:2) , , . . . , + √ m (cid:3) T ∈ R m . Theorem 2 (EDM characterization [19, Th. 3.3]).
A non-negative, symmetric,hollow matrix X ∈ R m × m is a Euclidean distance matrix if and only if the block (cid:98) X ∈ R ( m − ) × ( m − ) in Q ( − X ) Q = (cid:20) (cid:98) X dd T δ (cid:21) (1.5) is positive semi-definite. In this case, X is irreducibly embeddable in R q whereq = rank ( (cid:98) X ) ≤ m − . The problem of low-rank Euclidean distance matrix completion can now be for-mulated. Let D denote a partial Euclidean distance matrix, with entry D i j knownwhenever ( i , j ) ∈ Ω for some index set Ω , which is embeddable in R q . Without lossof generality, we make the following three simplifying assumptions on the partialmatrix D and index set Ω .1. (non-negative) D ≥ i.e., D i j ≥ i , j = , , . . . , m ); When two amino acids form a peptide bond, a water molecule is formed. An amino acid residue is what remains of each amino acid after this reaction. J.M. Borwein, M.K. Tam
2. (hollow) D ii = ( i , i ) ∈ Ω for i = , , . . . , m ;3. (symmetric) ( i , j ) ∈ Ω ⇐⇒ ( j , i ) ∈ Ω , and D i j = D ji for all ( i , j ) ∈ Ω .We define two constraint sets C = (cid:8) X ∈ S m : X ≥ , X i j = D i j for all ( i , j ) ∈ Ω (cid:9) , C = (cid:40) X ∈ S m : Q ( − X ) Q = (cid:20) (cid:98) X dd T δ (cid:21) , (cid:98) X ∈ S m − + , d ∈ R m − rank (cid:98) X ≤ q , δ ∈ R (cid:41) . (1.6)In light of Theorem 2, the problem of low-rank Euclidean distance matrix comple-tion can be cast as the two-set feasibility problemfind X ∈ C ∩ C . That is, a matrix X is a low-rank Euclidean distance matrix which completes D ifand only if X ∈ C ∩ C . Some comments regarding the constraint sets in (1.6) arein order.The set C encodes the experimental data obtained from NMR, and the a priori knowledge that the matrix is non-negative, symmetric and hollow. Its projection hasa simple formulae, as we now show. Proposition 1 (Projection onto C ). Let X ∈ R m × m . Then P C X is given element-wise by ( P C X ) i j = (cid:40) D i j , ( i , j ) ∈ Ω max { , X i j } , ( i , j ) (cid:54)∈ Ω for i , j = , , . . . , m . Proof.
Let Y be any matrix in C . We have (cid:107) X − Y (cid:107) F = ∑ ( i , j ) ∈ Ω ( X i j − Y i j ) + ∑ ( i , j ) (cid:54)∈ Ω s.t. X ij < ( X i j − Y i j ) + ∑ ( i , j ) (cid:54)∈ Ω s.t. X ij ≥ ( X i j − Y i j ) = ∑ ( i , j ) ∈ Ω ( X i j − D i j ) + ∑ ( i , j ) (cid:54)∈ Ω s.t. X ij < X i j + ∑ ( i , j ) (cid:54)∈ Ω s.t. X ij ≥ ( X i j − Y i j ) . (1.7)Let P be the matrix given by the proposed projection formula (clearly P ∈ C ). Then ∑ ( i , j ) (cid:54)∈ Ω s.t. X ij ≥ ( X i j − Y i j ) ≥ ∑ ( i , j ) (cid:54)∈ Ω s.t. X ij ≥ ( X i j − X i j ) = ∑ ( i , j ) (cid:54)∈ Ω s.t. X ij ≥ ( X i j − P i j ) . (1.8)By combining (1.7) and (1.8) we see that (cid:107) X − Y (cid:107) F ≥ (cid:107) X − P (cid:107) F for all Y ∈ C . Since C is closed and convex, P is the unique nearest point to X in C . (cid:117)(cid:116) Reflection methods for inverse problems 9
Remark 3.
Since C is a closed convex set, an alternative (less direct) proof of Propo-sition 1 can be given using the standard variational characterization of convex pro-jections [15, Th. 1.2.4].Using the necessary condition given by Theorem 2, the non-convex set C en-codes the a priori knowledge that the matrix of interest is a EDM together with thedimension of the space in which the corresponding points generating the matrix arecontained. We now derive the projection onto C . Theorem 3 (Nearest low-rank EDMs [3]).
Let X ∈ S m be a non-negative, hollowmatrix. ThenP C ( X ) = (cid:26) − Q (cid:20) (cid:98) Y dd T δ (cid:21) Q : Q ( − X ) Q = (cid:20) (cid:98) X dd T δ (cid:21) , (cid:98) X ∈ R ( m − ) × ( m − ) , d ∈ R m − , δ ∈ R , (cid:98) Y ∈ P M (cid:98) X (cid:27) , where M is the set of positive semi-definite matrices with rank q or less. In particular,P C ( X ) is a singleton if and only if P M (cid:98) X is a singleton.Proof.
Let Y be any matrix in C . That is, Y = (cid:20) (cid:98) Y cc T β (cid:21) , for some c ∈ R m − , β ∈ R , (cid:98) Y ∈ S . Using the orthogonality of Q , we compute (cid:107) X − Y (cid:107) F = (cid:107) Q ( X − Y ) Q (cid:107) F = (cid:107) Q ( − X ) Q − Q ( − Y ) Q (cid:107) F = (cid:13)(cid:13)(cid:13)(cid:13)(cid:20) (cid:98) X dd T δ (cid:21) − (cid:20) (cid:98) Y cc T β (cid:21)(cid:13)(cid:13)(cid:13)(cid:13) F = (cid:13)(cid:13)(cid:13)(cid:13)(cid:20) (cid:98) X − (cid:98) Y ( d − c )( d − c ) T ( δ − β ) (cid:21)(cid:13)(cid:13)(cid:13)(cid:13) F = (cid:107) (cid:98) X − (cid:98) Y (cid:107) F + (cid:107) d − c (cid:107) + | γ − β | . (1.9)To complete the proof we observe that (1.9) is minimized if and only if c = d , γ = β and (cid:98) Y ∈ P M (cid:98) X . (cid:117)(cid:116) The set M in Theorem 3 is a set of low-rank positive semi-definite matrices. Onemethod to compute its projection (and the one we will use) is by exploiting theeigen-decomposition of (cid:98) X . Denote by diag ( λ ) the diagonal matrix given by placingthe elements of the vector λ ∈ R m along the main diagonal. Let (cid:98) X = U diag ( λ ) U T be an eigen-decomposition (of (cid:98) X ) with λ ≥ λ ≥ · · · ≥ λ + q ≥ · · · ≥ λ m . A projection onto the set is then given by U diag (( λ + , λ + , . . . , λ + q , . . . , , )) U T , where x + denotes max { , x } . We apply the formulation of Section 1.6 to six proteins, shown in Table 1.1, obtainedfrom the
RCSB Protein Data Bank . As part of [3], reconstructions of the same sixproteins were attempted using a partial EDM containing only distances less than6 ˚A. Here we attempt reconstructions using partial EDMs which, in addition to theseshort-range distances, incorporate other a priori information. In particular, we in-clude inter-atomic distances greater than 6 ˚A for atoms from within the same residuein the partial EDM. This is reasonable since the structure of the individual residuesis known. For 1PTQ, this information gives approximately a further 0 .
2% of thetotal non-zero inter-atomic distances.
Table 1.1
Number of atoms, residues, known, and total non-zero inter-atomic distances in our sixtest proteins.Protein Atoms Residues Total Non-Zero Distances Known Non-Zero Distances1PTQ 404 50 81,406 8.9207%1HOE 581 74 168,490 6.4105%1LFB 641 99 205,120 5.6362%1PHT 988 85 236,328 4.6501%1POA 1067 118 568,711 3.6375%1AX8 1074 146 576,201 3.5606%
Our experiments were implemented in
Cython and performed on a machine hav-ing an Intel Xeon E5540 @ 2 . x , was converted to points z , z , . . . , z m ∈ R using Algorithm 1.4. Remark 4.
It is worth emphasing that our primary concern is the quality of the re-construction, rather than the time required to perform the reconstruction. This isbecause, if done well, one only needs to determine the conformation once.We report two error metrics, which we now explain. Denote the actual EDM by x actual . The first error metric is a measure of the error in the reconstructed EDM , andis given by EDM-error = (cid:107) x actual − x (cid:107) F = (cid:115) m ∑ i , j = (cid:12)(cid:12)(cid:12) x actual i j − x i j (cid:12)(cid:12)(cid:12) . RCSB Protein Data Bank:
Reflection methods for inverse problems 11
Fig. 1.4
Conversion of EDM to points in R q . Input : x ∈ X ; /* a Euclidean distance matrix */ L = I − ee T / n where e = ( , ,..., ) T ; τ = − LDL / USV T = SingularValueDecomposition ( τ ) ; Z = first q columns of U √ S ; z i = i th row of Z for i = , ,..., m ; Output : z , z ,..., z q ∈ R q ; /* points corresponding to x */ Denote the actual atom positions by z actual1 , z actual2 , . . . , z actual m ∈ R . The seconderror metric measures the error in the reconstructed atom positions z , z , . . . , z m ∈ R . Since EDMs are invariant under translation, reflection, and rotation of the pointsby which they are induced, we first perform a Procrustes analysis [16] to obtain (cid:101) z , (cid:101) z , . . . , (cid:101) z m ∈ R . These points are a best fit of the reconstructed points when theaforementioned transformations are allowed. The second error metric is given byPosition-error = (cid:115) m ∑ k = (cid:107) z actual k − (cid:101) z k (cid:107) . When comparing the relative size of these two errors, it is worth noting that the sum-mation in the EDM-error contains m terms whereas the summation in the position-error contains only 3 m . Remark 5 (Decibel error).
It is also common to consider the relative error in decibels(dB) , as was reported in [3]. That is,Relative error (dB) =
10 log (cid:18) (cid:107) P B R A x − P A x (cid:107) F (cid:107) P A x (cid:107) F (cid:19) . In this study the relative error in decibels is not reported. This is unnecessary be-cause the stopping criterion used in Algorithm 1.3 is equivalent to requiring thatthe decibel error be less than 10 log ( ε ) . Requiring that ε = − corresponds toaiming at a relative error of − Remark 6 (Stopping criterion and tolerance).
In the computational experiments thatfollow, the stopping tolerance is taken to be ε = − . We now provide some justi-fication for this choice.For each of the six proteins, Figure 1.5 shows the relative error as a function ofthe number of iterations starting from a given initial point for the Douglas–Rachfordmethod. • When the number of iteration is less than 5000 the relative error exhibits non-monotone oscillatory behaviour — which seems to provide much of the potency
Fig. 1.5
The relative error as a function of iterations (vertical axis is logarithmic). of the method. It seems to allow the reflection method to sample regions andavoid settling at an inferior local minimum of the configuration space. In [3]we observed that the alternating projection method, which is monotonic, fails toproduce good reconstructions. • When the relative error is between 10 − and 10 − , it decreases sharply afterwhich a period of more predictable decrease is observed. • Beyond this point slower progress is made. We therefore choose our stoppingtolerance to be ε = − so that the algorithm will terminate in this region.The change in successive iterates was found to also exhibit similar behavior (notshown), so is another suitable candidate for a stopping criterion.It is worth noting that there are many other techniques for solving (variants of) theprotein conformation problem (see for instance [21]). Such a discussion, however,is beyond the scope of this chapter. Table 1.2 gives results for the basic Douglas–Rachford algorithm presented in Al-gorithm 1.3. We make some comments regarding these results.The EDM-error increases with increasing problem size; yet the same trend isnot observed for the position-error for which 1PHT reported the largest error. Forall of the proteins studied, the differences between the average and worst case re-sults for the position-errors were small. This strongly suggests that the method canconsistently produce a EDM which gives the desired atomic positions.
Reflection methods for inverse problems 13
The second column of Figure 1.6 shows the conformation of the basic Douglas–Rachford reconstructions, which are visually indistinguishable from the actual con-formation shown in the first column. This is an improvement from what was reportedin [3] whose Douglas–Rachford reconstructions of two of the larger proteins, 1POAand 1AX8, gave unrealistic conformations consisting of disjoint blocks of atoms. Inlight of Remark 6 it is likely that this was due to premature algorithm termination.
Table 1.2
Average (worst) results from five random replications of the basic Douglas–Rachfordalgorithm with ε = − .Protein EDM-Error Position-Error Iterations Time (h)1PTQ 3 . ( . ) . ( . ) . ( ) . ( . ) . ( . ) . ( . ) . ( ) . ( . ) . ( . ) . ( . ) . ( ) . ( . ) . ( . ) . ( . ) . ( ) . ( . ) . ( . ) . ( . ) . ( ) . ( . ) . ( . ) . ( . ) . ( ) . ( . ) In our formulation of the protein confirmation problem, the most expensive stepis the computation of the projection onto the rank constraint C . Thus requires theeigen-decomposition of a ( m − ) × ( m − ) symmetric matrix. In this section wepropose problem specific heuristics which allow for this computation to sometimesbe avoided.One idea to avoid performing the eigen-decomposition is to not update the se-quence ( r n ) ∞ n = in Algorithm 1.3 at every iteration but only periodically. This ap-proach is described in Algorithm 3, and results, with updates only every third time,in Table 1.3.We now compare the results of this section to those of Section 1.7.1. A small in-crease in the position-errors, and a larger increase in the EDM-errors was observed.The number of iterations required also increased, with this number almost doublingfor 1PTQ. For all six test proteins, the total time required was less. The biggestimprovement was 1POA whose total time was more than halved. The quality of thereconstructed conformations seem not to be adversely effected by the use of periodicrank projections, as can be seen in Figure 1.6. P T Q HO E L F B P H T P OA AX Fig. 1.6
The conformations of the six proteins, and their three Douglas–Rachford reconstructions.
Fig. 1.7
The Douglas–Rachford algorithm with T -periodic projections onto the set B . Input : x ∈ X , T ∈ N and ε > n = p ∈ P A ( x ) ; while n = or (cid:107) r n − p n (cid:107) > ε (cid:107) p n (cid:107) doif n mod T = then r n ∈ P B ( p n − x n ) ; else r n = r n − ; end x n + = x n + r n − p n ; p n + ∈ P A ( x n + ) ; n = n + endOutput : p n ∈ X Reflection methods for inverse problems 15
Table 1.3
Average (worst) results from five random replications of the Douglas–Rachford algo-rithm with periodic rank projections with T = ε = − .Protein EDM-Error Position-Error Iterations Time (h)1PTQ 4 . ( . ) . ( . ) . ( ) . ( . ) . ( . ) . ( . ) . ( ) . ( . ) . ( . ) . ( . ) . ( ) . ( . ) . ( . ) . ( . ) . ( ) . ( . ) . ( . ) . ( . ) . ( ) . ( . ) . ( . ) . ( . ) . ( ) . ( . ) In Sections 1.7.1 & 1.7.2 we considered the physically realistic setting in whichdistances below the threshold of 6 ˚A were known. As noted, when the number ofatoms in a protein increases, the proportion of inter-atomic distances below thisthreshold compared to the total number of (non-zero) distances decreases.To better understand the Douglas–Rachford method applied to larger probleminstances, we performed the same reconstruction as in Section 1.7.1 but with thepercentage of known non-zero distances constant. More precisely, we assumed thatthe smallest 10% of inter-atomic distances were known.
Table 1.4
Average (worst) results from five random replications of the basic Douglas–Rachfordalgorithm from the smallest 10% of inter-atomic distances with ε = − .Protein EDM-Error Position-Error Iterations Time (h)1PTQ 3 . ( . ) . ( . ) . ( ) . ( . ) . ( . ) . ( . ) . ( ) . ( . ) . ( . ) . ( . ) . ( ) . ( . ) . ( . ) . ( . ) . ( ) . ( . ) . ( . ) . ( . ) . ( ) . ( . ) . ( . ) . ( . ) . ( ) . ( . ) As could perhaps be predicted, when more distance information is incorporatedthe error metrics, and the number of iterations decrease. Problem size and EDM-error do not correlate as strongly compared to the results of Section 1.7.1. However,the general trend that larger problem sizes give larger EDM-errors is still observed.The most notable improvement, when compared to Section 1.7.1, is the position-error for 1PHT. This suggests that in the realistic setting of Section 1.7.1 the un-derlying protein’s conformation ( e.g., a compact or a dispersed conformation) is animportant factor in the difficulty of the reconstruction problem.
Ionic liquids (ILs) are salts ( i.e., they are comprised of positively and negativelycharged ions) having low melting points, typically occupying the liquid state at roomtemperature. An analogous reconstruction problem arising in the context of ionicliquid chemistry is to determine a given ionic liquid’s bulk structure . That is, theconfiguration of its ions with respect to each other (the structure of the individualions is known).In this section, we applied the Douglas–Rachford method to a simplified versionof this problem. Entries of the partial EDM are assumed to be known whenever thetwo atoms are bonded ( i.e., when their
Van der Waals radii taken from [8] overlap).Table 1.5 reports results for a propylammonium nitrate (PAN) data set consistingof 180 atoms. The corresponding rank-3 EDM completion problem has a total of32,220 non-zero inter-atomic distances of which 5.95% form the partial EDM.
Table 1.5
Average (worst) results from five random replications of the basic Douglas–Rachfordalgorithm, applied to ionic liquid bulk structure determination, with ε = − .EDM-Error Position-Error Iterations Time (h)0 . ( . ) . ( . ) . ( ) . ( . ) As was the case in the protein conformation application, the difference betweenthe average and worst case results for the two error metrics is observed to besmall. The actual conformation of PAN, and its Douglas–Rachford reconstructionare shown in Figure 1.8. A high degree of visual coincidence is observed, althougha small amount of the finer detail is missing.
Fig. 1.8
The actual conformation (left) and Douglas–Rachford reconstruction (right) of PAN. Notethe two poorly reconstructed hydrogen atoms (white) in the left configuration. Reflection methods for inverse problems 17
We have shown that the Douglas–Rachford reflection method can successfullysolve the protein conformation determination problem by directly addressing a non-convex matrix completion problem. This is also the case for an analogous ionic liq-uid bulk structure determination problem. It is worth emphasising again that the cur-rent literature provides no theoretical justification for the method to work at all, letalone so well. Modifications of the method have also been shown to reduce computa-tional times without significantly effecting the quality of the results. This promisingdemonstration of the method begs further attention, both in improving theoreticalunderstanding, and in the refinement and investigation of these and further applica-tions.
Acknowledgements.
The authors wish to thank Dr Alister Page for introducing usto the bulk structure determination problem, and for kindly sharing the PAN dataset. The work of JMB is supported, in part, by the Australian Research Council.The work of MKT is supported, in part, by an Australian Postgraduate Award.
References
1. Arag´on Artacho, F., Borwein, J.: Global convergence of a non-convex Douglas–Rachford it-eration. J. Glob. Optim. (3), 753–769 (2013).2. Arag´on Artacho, F., Borwein, J., Tam, M.: Recent results on Douglas–Rachford methods forcombinatorial optimization problems. J. Optim. Theory Appl. (in press, 2013).3. Arag´on Artacho, F., Borwein, J., Tam, M.: Douglas–Rachford feasibility methods for matrixcompletion problems. ANZIAM J. (in press, 2014).4. Bauschke, H., Bello Cruz, J., Nghia, T., Phan, H., Wang, X.: The rate of linear convergenceof the Douglas–Rachford algorithm for subspaces is the cosine of the Friedrichs angle. J.Approx. Theory 185, 63–79 (2014).5. Bauschke, H., and Combettes, P.: Convex analysis and monotone operator theory in Hilbertspace. Springer, New York (2011).6. Bauschke, H., Combettes, P., Luke, D.: Finding best approximation pairs relative to two closedconvex sets in Hilbert spaces. J. Approx. Theory (2), 178–192 (2004).7. Bauschke, H., Noll, D., Phan, H.: Linear and strong convergence of algorithms involving av-eraged nonexpansive operators. arXiv preprint arXiv:1402.5460 (2014).8. Bondi, A.: Van der Waals Volumes and Radii. J. Phys. Chem. (3):441–51 (1964).9. Borwein, J., Lewis, A.: Convex analysis and nonlinear optimization. Springer (2006).10. Borwein, J., Sims, B.: The Douglas–Rachford algorithm in the absence of convexity. In: Fixed-Point Algorithms for Inverse Problems in Science and Engineering, pp. 93–109. Springer(2011).11. Borwein, J., Tam, M.: The cyclic Douglas–Rachford method for inconsistent feasibility prob-lems. J. Nonlinear Convex Analysis, accepted March 2014. arXiv preprint arXiv:1310.2195(2013).12. Borwein, J., Tam, M.: A cyclic Douglas–Rachford iteration scheme. J. Optim. Theory Appl. (1), 1–29 (2014).13. Borwein, J., Zhu, Q.: Techniques of Variational Analysis, CMS Books in Mathematics , vol. 20.Springer-Verlag, New York (2005, Paperback, 2010).14. Berman, A., Shaked-Monderer, N.: Completely positive matrices. World Scientific, Singapore(2003).8 J.M. Borwein, M.K. Tam15. Cegielski, A.: Iterative methods for fixed point problems in Hilbert space,
Lecture Notes inMathematics , vol. 2057. Springer, London (2012).16. Dattorro, J.: Convex optimization & Euclidean distance geometry. Meboo Publishing USA(2005).17. Elser, V., Rankenburg, I., Thibault, P.: Searching with iterated maps. Proc. Natl. Acad. Sci. (2), 418–423 (2007).18. Gravel, S., Elser, V.: Divide and concur: A general approach to constraint satisfaction. Phys.Rev. E (3), 036,706 (2008).19. Hayden, T., Wells, J.: Approximation by matrices positive semidefinite on a subspace. LinearAlgebra Appl. , 115–130 (1988).20. Hesse, R., Luke, D.: Nonconvex notions of regularity and convergence of fundamental algo-rithms for feasibility problems. SIAM J. Optim. (4), 2397–2419 (2013).21. Seo, J., Kim, J.-K., Ryu, J., Lavor, C., Mucherino, A., and Kim, D.-S.: BetaMDGP: Proteinstructure determination algorithm based on the Beta-complex. Trans. Comput. Sc.8360