Completely positive factorizations associated with Euclidean distance matrices corresponding to an arithmetic progression
Damjana Kokol Bukovšek Thomas Laffey Helena Šmigoc































lab1,lab2]Damjana Kokol Bukovšekfootnote
[footnoste]Damjana Kokol Bukovšek acknowledges financial support from the Slovenian Research Agency (research core funding No. P1-0222). The work on the paper was done during her three-month visit to University College Dublin.
lab3]Thomas Laffey
lab3]Helena Šmigoc
[lab1]School of Economics and Business, University of Ljubljana, Slovenia
[lab2]Institute of Mathematics Physics and Mechanics, Ljubljana, Slovenia
[lab3]School of Mathematics and Statistics, University College Dublin, Ireland


Euclidean distance matrices corresponding to an arithmetic progression have rich spectral and structural properties. We exploit those properties to develop completely positive factorizations of translations of those matrices. We show that the minimal translation that makes such a matrix positive semidefinite results in a completely positive matrix. We also discuss completely positive factorizations of such matrices over the integers. Methods developed in the paper can be used to find completely positive factorizations of other matrices with similar properties. 














completely positive matrices Euclidean distance matrices 


15B48 15A23 15B36





























































 


 
















§  INTRODUCTION 


An n× n real symmetric matrix A is completely positive, if it can be written as A=BB^T for some n × k entry-wise non-negative matrix B. The minimum number of columns k in such a non-negative factor B is called the cp-rank of A. The set of completely positive matrices forms a closed, convex cone, and better understanding of this cone is not only an interesting question, but it also plays an important role in the field of copositive optimization, <cit.>. 

It is immediate from the definition, that any completely positive matrix is non-negative and positive semidefinite. A matrix that is both non-negative and positive semidefinite is said to be doubly non-negative. Not every n × n doubly non-negative matrix of order n ≥ 5 is completely positive, and deciding if a given matrix is completely positive is a hard problem that has been a topic of intense research over the years. For a comprehensive survey on completely matrices we refer to <cit.>. 

Research has shown that the question can be solved for different classes of matrices. Here we mention two representative results that are different in flavour. In <cit.> it is shown that every diagonally dominant non-negative matrix is completely positive. Kogan and Berman <cit.> introduced a graph theoretic view on the problem. For a given n× n symmetric matrix A, let G(A) be a graph on n vertices with an edge (i,j) if and only if a_ij≠ 0.  A graph G is said to be completely positive if every doubly non-negative matrix with G(A)=G is completely positive. In <cit.> it is shown that a graph G is completely positive if and only if it does not contain an odd cycle of length greater than 4.

Finding a factorization of a given completely positive matrix is a natural question in this context. A factorization algorithm that can factorize any matrix in the interior of the cone of completely positive matrices is offered in <cit.>, explicit factorization of diagonally dominant matrices is presented in <cit.>, and generalized to matrices whose comparison matrix is positive definite in <cit.>. Factorization of matrices with conditions on their graph is studied for example in <cit.>. Those factorizations are typically not optimal, the number of columns in the factorization matrix is larger than the minimal possible. 
Finding and optimizing the cp-rank adds additional complexity to an already hard problem. 

Once we know that a factorization of a given matrix A exists, we may want to impose further properties on the factors. A matrix A has a rational cp-factorization, if it can be decomposed as BB^T where B the entries of B are non-negative and rational. Integer cp-factorization is defined in a similar way.  Every rational matrix which lies in the interior of the cone of completely positive matrices has a rational cp-factorization <cit.>. Integer question seems to be harder to study, and has only recently been answered for 2 × 2 matrices <cit.>. For a detailed state-of-the art account on those two questions we refer the reader to <cit.>.

In this paper we study completely positive factorizations of translations of Euclidean distance matrices corresponding to an arithmetic progression. Distance matrices have rich structural and spectral properties, and have been extensively studied in the literature <cit.>. Our work was particularly motivated by results on factorizations of distance matrices <cit.>. In this paper we develop further factorizations of translations of distance matrices  corresponding to an arithmetic progression, that depend on specific properties of this class. 
 Methods, that we introduce, can be adapted to other matrices with similar properties. While Euclidean distance matrices are well behaved in many ways, they are not positive definite. The first factorization that  we develop has a factor with one negative element, and subsequent factorizations are developed for translations of  Euclidean distance matrices. In Section <ref> we find the best possible result for this question. Finally, in Section <ref> we present some factorizations, where the factors have integer entries. 

We will use standard notation in the paper. The n × n identity matrix will be denoted by I_n, and the n × n nilpotent Jordan block by J_n. B(i) will denote the matrix obtained from B by deleting the i-th row and column. More generally, for ℐ⊆{1,2,…,n}, we will denote by B(ℐ) the matrix obtained from B by deleting the row and columns indexed by ℐ. 



§ EUCLIDEAN DISTANCE MATRICES CORRESPONDING AN ARITHMETIC PROGRESSION


For a set of distinct real numbers 𝒮 = {a_1, a_2, ..., a_n}, we  define the Euclidean distance matrix to be the matrix EDM(𝒮)=(a_ij)_i,j=1^n with a_ij = (a_j - a_i)^2. Such a matrix EDM(𝒮) is clearly entrywise non-negative, and it has rank 3, see <cit.>.


In this work we are interested in a family of Euclidean distance matrices that correspond to an arithmetic progression, i.e. 𝒮 = {a, a+d, ..., a+(n-1) d} for some positive real numbers a and d, and in particular for  𝒮 = {1, 2, ..., n}. We will denote 
    A_n := EDM({1, 2, ..., n}),
 hence (A_n)_ij = (j-i)^2.
We can limit our discussion to this special case, as for 𝒮 = {a, a+d, ..., a+(n-1) d} we have EDM( S) = d^2 A_n, and all the properties of A_n are easily adapted to this more general case. In our first result we explicitly compute the eigenvalues and the eigenvectors of A_n. 


The nonzero eigenvalues of A_n are equal to:

    λ_1,2 = 1/12n(n^2-1) ±√(1/240n^2(n^2-1)(3n^2-7)), 
    λ_3 = -1/6n(n^2-1),

λ_1 > 0 > λ_2 > λ_3. The eigenvector corresponding to the eigenvalue λ_3 is equal to 
w = (w_i)_i=1^n with w_i:=n+1-2i for i=1,…, n. For j=1,…, n-3 let v(j) denote the vector with only four nonzero elements, defined as follows: v(j)_j=1, v(j)_j+1=-3, v(j)_j+2=3, v(j)_j+3=-1. Then the vectors v(j), j=1,…,n-3, form a basis of the null space of A_n.  



We can prove that w is an eigenvector for A_n corresponding to λ_3 by direct computation:

    (A_nw)_i = ∑_j=1^n (j-i)^2(n+1-2j) = -1/6n(n^2-1)(n+1-2i).


It is also straightforward to check that A_nv(j)=0, for j=1,…,n-3. 
Since rank of A_n is 3, we still need to compute two more nonzero eigenvalues. 





To compute the remaining two eigenvalues, we note that the trace of A_n is equal to zero, i.e. λ_1+λ_2+λ_3=0, 
and we compute the trace of A_n^2: 

    λ_1^2 + λ_2^2 + λ_3^2 = ∑_i=1^n ∑_j=1^n (j-i)^4 = 1/30n^2(n^2-1)(2n^2-3).

Taking into account that λ_3 is known, we obtain the following quadratic equation for λ_1 and λ_2:

    180 λ^2 - 30 n(n^2-1) λ - n^2(n^2-1)(n^2-4) = 0 .

Solving this equation gives us: 
    λ_1 = 1/12n(n^2-1) + √(1/240n^2(n^2-1)(3n^2-7))

and 
    λ_2 = 1/12n(n^2-1) - √(1/240n^2(n^2-1)(3n^2-7)).

Finally, the inequality λ_1>λ_2 > λ_3 is straightforward to check.


Since A_n has two negative eigenvalues it is clearly not completely positive. Still, the aim of this work is to develop factorizations of matrices A_n and A_n+α I_n that take into account not only the nonnegativity of A_n, but also as many of the following properties of A_n as possible: symmetry, integer entries, pattern, and low rank of A_n. We won't be able to attend to all these properties with a single factorization, and we will allow translations of A_n to make it positive semidefinite.  



§ STRAIGHTFORWARD FACTORIZATIONS RELATED TO A_N



First we offer a factorization of A_n, that is reminiscent of a completely positive factorization, but has an additional factor with a negative element. 


We can write A_n=L_nRL_n^T, where L_n is a non-negative n × 3 matrix and

    R =  0     1     1 
     1     -6     1 
     1     1     0 .




Notice that the vectors v(j), j=1,…,n-3, defined in Lemma <ref> are the first n-3 columns of the matrix (I_n-J_n^T)^3.
 It follows that first n-3 columns of the matrix 
A_n(I_n-J_n^T)^3 are equal to zero, and, by symmetry, the first n-3 rows of Â_n:=(I_n-J_n)^3A_n(I_n-J_n^T)^3 are also zero. 
Hence,  only the lower-right 3× 3 corner of Â_n is non-zero, and it is straightforward to check that this corner is equal to 

    R =  0     1     1 
     1     -6     1 
     1     1     0 .







It follows that A_n=(I_n-J_n)^-3(0_n-3⊕ R)(I_n-J_n^T)^-3, and L_n is equal to the last three columns of (I_n-J_n)^-3. The claim is proved by noting that (I_n-J_n)^-1 = I_n + J_n + J_n^2 + ... + J_n^n-1 is a non-negative matrix. 


An interesting feature of the factorization in Theorem <ref> is that the matrix R is independent of n. Since we cannot hope to find a completely positive factorization of A_n, we focus on the factorizations of matrices of the form 
    A_n+g(n) I_n.
 Clearly, A_n+g(n) I_n will become diagonally dominant, and hence completely positive, for all large enough g(n). The minimal g(n) for which A_n is diagonally dominant is easily computed, and is equal to
g_D(n):=∑_j=1^n-1 j^2=1/6n(n-1)(2n-1). On the other hand, the minimal g(n) for which A_n+g(n) I_n will become positive semidefinite can be deduced from Lemma <ref>:  f(n)= 1/6n(n^2-1). The minimal g(n) for which A_n+g(n)I_n is completely positive therefore satisfies the following inequality: 

    1/6n(n-1)(n+1) ≤ g(n) ≤1/6n(n-1)(2n-1).

The main result of this paper will show that the lower bound is achieved, but first we will improve the upper bound in (<ref>) using an inductive approach. 
Informed by (<ref>) we denote f(n)=1/6n(n-1)(n+1), and look for functions g(n) of the form g(n)=q f(n), for which we can prove that A_n+g(n)I_n is completely positive using a straightforward inductive proof. 

Note that A_n-1=A_n(1)=A_n(n), and the most natural inductive approach to the problem is to write   

    A_n + g(n)I_n=(A_n-1+g(n-1)I_n-1) ⊕ 0_1 +R_n,

where v_n-1=(n-1)^2     (n-2)^2    …    1^T and 
    R_n=(g(n)-g(n-1)) I_n-1    v_n-1
     v_n-1^T     g(n).

If we already know that A_n-1+g(n-1)I_n-1 is completely positive by inductive hypothesis, then it is sufficient to prove that R_n is completely positive to determine the complete positivity of A_n. Since the graph of R_n contains no odd cycles of order greater than 3, the matrix R_n is completely positive as soon as it is doubly non-negative. Note that to determine that R_n is positive semidefinite we only need to show that R_n >0 from consecutive minors result, see for example <cit.>. Inserting g(n)=q f(n) into 
 
    R_n=(g(n)-g(n-1))^n-2(g(n)(g(n)-g(n-1))-v_n-1^Tv_n-1),

 and collecting (g(n)(g(n)-g(n-1))-v_n-1^Tv_n-1) in terms of n, gives us q= 2
   √(3/5) as the optimal choice for q for this approach. 
  
Next we consider an induction step that takes us from n-2 to n. Note that A_n-2=A_n({1,2})=A_n({n-1,n})=A_n({1,n}), and it turns out that developing our induction process from A_n-2=A_n({1,n}) gives the best bound. 
Hence, we write 
 
    A_n + g(n)I=0_1 ⊕( A_n-2+g(n-2)I_n-2) ⊕ 0_1 + Q_n,

where 

    Q_n = g(n)     u_n-2^T     (n-1)^2 
     u_n-2    (g(n)-g(n-2))I_n-2    v_n-2
     (n-1)^2     v_n-2^T     g(n) ,

u_n-2 =  1     2^2     ...     (n-2)^2^T, and  v_n-2 = (n-2)^2     ...     2^2     1^T.
As above, it is sufficient to consider the complete positivity of Q_n. The graph of Q_n contains no odd cycles of order greater than 3, so Q_n will be completely positive as soon as it will be 
positive semidefinite. To study the positive definiteness of Q_n, we consider a matrix that is permutationally similar to Q_n:
 
    Q̂_n= (g(n)-g(n-2))I_n-2    u_n-2    v_n-2
     
                                           u_n-2^T     g(n)     (n-1)^2 
    
                                           v_n-2^T     (n-1)^2     g(n).

Assuming g(n) >g(n-2), Q̂_n will be positive semidefinite if and only if Q̂_n(n) >0 and Q̂_n≥ 0. Using Schur complement to compute those determinants, and considering the leading coefficient in n, we determine that Q̂_n(n) >0 for all q > √(6/5), and that Q̂_n>0 if and only if 
 
    δ_n g(n)-u_n-2^Tu_n-2   δ_n(n-1)^2-v_n-2^Tu_n-2
    δ_n(n-1)^2-v_n-2^Tu_n-2   δ_n g(n)-v_n-2^Tv_n-2>0,

where δ_n=(g(n)-g(n-2)). Taking into account:

    u_n-2^Tu_n-2 =v_n-2^Tv_n-2   =1/30  (n-1)(n-2) (2 n-3) (3 n^2-9 n+5),  
    
    v_n-2^Tu_n-2   =1/30 n(n-1)(n-2) (n^2-2 n+2 ),

and g(n)=q f(n), we determine that q= √(7/5) is optimal for this method. This gives us the following proposition that strengthens the bound given in  (<ref>). 


The matrix A_n+√(7/5) f(n)I_n is completely positive. 




§  OPTIMAL G(N) 


At this point we would like to strengthen Proposition <ref> to q=1. To this end we denote B_n=A_n+f(n)I_n. Note that B_n=0, hence we are operating at the boundary of positive semidefinite matrices. The following simple and well known result points out an additional property that the factors in the cp factorization of a semidefinite matrix must satisfy. 


Let A=BB^T be a completely positive matrix with a nontrivial kernel 𝒩. Then B^Tv=0 for every v ∈𝒩.
 

The proof of the main result in this section will be rather technical, as we will explicitly develop a completely positive factorization of B_n. We will aim to write B_n as the sum of rank one completely positive matrices, and a matrix with particular structure that is made explicit in the following lemma. 


Let D_1 ∈ M_m() and D_2 ∈ M_k() be diagonal matrices, and let C be an m × k non-negative matrix. Then any matrix of the form

    A =  D_1     C 
     C^T     D_2∈ M_m+k(),

with an eigenvector w=w_1^T     -w_2^T^T corresponding to a non-negative eigenvalue, where w_1 ∈ M_m() and w_2 ∈ M_k() are strictly positive, is completely positive. 








First we observe that the matrix B=(I_m ⊕ -I_k)A(I_m ⊕ -I_k) has nonpositive off-diagonal elements, and a positive eigenvector ŵ=w_1^T     w_2^T^T corresponding to a non-negative eigenvalue. We deduce that Bŵ≥ 0, and hence B is an M-matrix, see for example <cit.>. Since B is also symmetric, it is positive semidefinite. This, in particular, proves that D_1 and D_2 are non-negative. Since A is similar to B, we now know that A is doubly non-negative. Moreover, the graph of A is bipartite, hence A is completely positive, as it is proved in <cit.>. 
















Proposition <ref> and Lemma <ref> bring to light, why it is convenient to know the singular vector of B_n explicitly when we are developing the completely positive factorization of B_n. 
In what follows let us denote the eigenvector corresponding to the eigenvalue λ_3 given in Lemma <ref> as w(n), or just w, when n is clear from the context. (Clearly, w(n) is a singular vector for B_n.) Part of the proof will be split into an even and an odd case. One of the reasons for this is the slight difference in the pattern of w(n) in those cases as w(2m)=v_2m,-v_2m'^T and w(2m+1)=v_2m+1,0,v_2m+1'^T, where v_n=n-1     n-3    …    n-(2 ⌊n/2⌋-1) ^T and v_n'= n-(2 ⌊n/2⌋-1)   …    n-3     n-1 ^T.



Let K_n ∈ M_n() be the matrix with ones on the anti-diagonal and zeros elsewhere. Then K_n^2=I_n, K_nA_nK_n=A_n and K_nv_n=v_n'. 



Let n=2m. There exists a non-negative matrix U satisfying U^Tw(n)=0, diagonal matrices D, D' ∈ M_m() and a non-negative matrix C ∈ M_m(), such that

    B_n = UU^T + D     C 
     C^T     D'.



 
As suggested by the result we will view B_n as a 2 × 2 block matrix: 

    B_n =  A_m + f(n)I_m     H_m 
     H_m^T     A_m + f(n)I_m,

where H_m=(h_ij) ∈ M_m() is defined by
h_ij=(m+j-i)^2. We will construct U by finding two matrices V and V' satisfying 
 UU^T=VV^T+V'V'^T,
 where each column of V corresponds to an off-diagonal pair of entries in the upper-left block of B_n and each column of V' corresponds to an off-diagonal pair of entries in the lower-right block of B_n. 

Explicitly, V:=∑_i=1^m-1∑_j=i+1^m u_ij u_ij^T,
where 

    u_ij   := (j-i) (e_i + e_j + α_ij e_m+k),  
     k:= ⌊m+j/2⌋ + 1 - i,  and 
    α_ij   := 2/2k-1(2m-i-j+1).

Note that the matrix u_ij u_ij^T has exactly nine nonzero entries, three on the main diagonal, one at the position (i,j), two in the upper-right block and symmetrically three more below the main diagonal. The entry in (i,j) position is equal to (i-j)^2, and hence equal to the (i,j)-th entry of B_n. Furthermore, once k is chosen, the value of α_ij is determined from the condition u_ij^Tw(n)=0. (While there is some choice in what k can be, we need to make sure that we choose it in such a way that the off-diagonal blocks remain non-negative after subtraction from B_n.) Thus the matrix VV^T has the form

    VV^T = A_m + E_1     S 
     S^T     E_2 ,

where E_1 and E_2 are diagonal matrices, and satisfies the condition V^Tw(n)=0. 

We define V':=K_nV. Using block partition 
    K_n=0_m     K_m 
     K_m     0_m
 and Lemma <ref> it is easy to determine that

    V'V'^T= E_1'     S' 
     S'^T     A_m+E_2'  and  V'^Tw(n)=0,

where E_1'=K_mE_2K_m, E_2'=K_mE_1K_m and S'=K_mS^TK_m. In particular, 
S'_xy=S_m+1-y,m+1-x.













To complete the proof, we need to show that the matrix C=H_m-S-S' is non-negative. 
 Let us first consider the (x, y)-entry in S: s_xy=∑_1 ≤ i < j ≤ m(u_iju_ij^T)_x,m+y.
 The matrix u_iju_ij^T has a nonzero entry at (x, m+y) if and only if (x, y) = (i, k) or (x, y) = (j, k).
In the first case we have i=x and y = ⌊m+j/2⌋ + 1 - i, thus j= 2(x+y-1)-m or j= 2(x+y-1)+1-m.
In the second case we have j=x and i = ⌊m+x/2⌋ + 1 - y. So, at most three 
columns of matrix VV^T contribute to the entry s_xy. We will refer to this contributions as matrices S_I, S_II and S_III, whose nonzero entries are defined as follows: 

    (S_I)_xy   =(u_x,2(x+y-1)-mu_x,2(x+y-1)-m^T)_x,y+m for 1 ≤ x < 2(x+y-1)-m ≤ m 
    
    (S_II)_xy   =(u_x,2(x+y-1)+1-mu_x,2(x+y-1)+1-m^T)_x,y+m for 1 ≤ x < 2(x+y-1)+1-m ≤ m 
    
    (S_III)_xy   =(u_⌊m+x/2⌋+1-y,xu_⌊m+x/2⌋+1-y,x^T)_x,y+m for 1 ≤⌊m+x/2⌋+1-y < x ≤ m.










The computation of relevant entries of u_iju_ij^T and reordering of the bounds, gives us the following formulas for the nonzero entries of S_I, S_II and S_III:

    (S_I)_xy   = 2/2y-1(x+2y-m-2)^2(3m-3x-2y+3)  for   m+3-2y ≤ x ≤ m+1-y, 
    
    (S_II)_xy   = 2/2y-1(x+2y-m-2)^2(3m-3x-2y+2)  for   m+2-2y ≤ x ≤ m-y, 
     
    (S_III)_xy   =2/2y-1(x-⌊m+x/2⌋+y-1)^2(2m-⌊m+x/2⌋+y-x) 
        for max{m+3-2y,2y-m-1}≤ x ≤ m.

Note that outside the indicated regions the entries of those matrices are equal to zero.










To simplify the calculations we estimate the matrix S_III by a matrix Ŝ_III with nonzero entries in positions (x,y) satisfying m+3-2y ≤ x ≤ m and defined by: 

    (S_III)_xy≤ (Ŝ_III)_xy = 2/2y-1(x-m+x-1/2+y-1)^2(2m-m+x-1/2+y-x).


Since the matrices H_m and S + S' are symmetric with respect to the counter-diagonal it is enough to prove the nonnegativity of elements (H_m-S-S')_xy that satisfy x+y ≤ m+1. Inside this region we will consider four subregions, that are determined by the location of nonzero entries of matrices S_I, S_II, S_III, S'_I, S'_II and S'_III. Before we define the regions, we observe that (S_II)_xy=0 for x+y≥ m+1, hence (S'_II)_xy=0 for x+y≤ m+1 and S'_II doesn't need to be considered. Similarly, (S_I)_xy=0 for x+y ≥ m+2, so (S'_I)_xy only contribution is when x+y=m+1. At this point we also note that (S)_m1=(S')_m1=0. Now let us define and study each of the four regions:
 


  * {(x,y);  x ≤ m+1-2y}, only (Ŝ'_III)_xy≠ 0, hence C_xy=(H_m-S'_III)_xy≥ (H_m-Ŝ'_III)_xy.  Introducing new variables u=y-1≥ 0 and v=m+1-x-2y≥0 gives us:

    C_xy≥3 + 55 u + 129 u^2 + 81 u^3 + 2 v + 52 u v + 66 u^2 v + 12 u v^2/4 (3 + 4 u + 2 v) > 0.



  * {(x,y);  x = m+2-2y, y ≥ 2)}, (S_II)_xy≠ 0 and (Ŝ'_III)_xy≠ 0, hence C_xy≥ (H_m-S_II-Ŝ'_III)_xy, which we express in terms of u=y-2≥ 0. We get

    C_xy≥320 + 1184 u + 1558 u^2 + 855 u^3 + 162 u^4/4 (3 + 2 u) (5 + 4 u) > 0.



  * {(x,y);  m+3-2y ≤ x ≤ m-y }, (S_I)_xy≠ 0, (S_II)_xy≠ 0, (Ŝ_III)_xy≠ 0, (Ŝ'_III)_xy≠ 0, and C_xy≥ (H_m-S_I-S_II-Ŝ_III-Ŝ'_III)_xy.
 Change of variables u=x+2y-m-3≥ 0 and v=m-x-y≥0 gives us

    C_xy≥p(u,v)/4 (5 + 2 u + 2 v) (7 + 2 u + 4 v) > 0,

where

    p(u,v)    = 1261 + 1487 u + 601 u^2 + 122 u^3 + 12 u^4 + 3629 v + 3306 u v  + 
       +  849 u^2 v + 72 u^3 v + 3565 v^2 + 2351 u v^2 + 312 u^2 v^2 + 1371 v^3  + 
       +  516 u v^3 + 162 v^4.



  * {(x,y);   x=m+1-y, y ≥ 2 }, (S_I)_xy=(S'_I)_xy≠ 0, (Ŝ_III)_xy=(Ŝ'_III)_xy≠ 0, and C_xy≥ (H_m-2S_I-2 Ŝ_III)_xy, 
which we express in terms of u=y-2≥ 0 to get

    C_xy≥6 + 16 u + 12 u^2 + 3 u^3/2 (3 + 2 u) > 0.







































Let n=2m+1. There exists a non-negative matrix U satisfying U^Tw(n)=0, α≥ 0, diagonal matrices 
D, D' ∈ M_m() and a non-negative matrix C ∈ M_m(), such that

    B_n = UU^T + D     0     C 
     0    α    0 
     C^T     0     D'.




















The proof will go along the same lines as the proof of Lemma <ref>. This time we will view B_n as a 3 × 3 block matrix: 

    B_n =  A_m + f(n)I_m     b_m     H_m 
     b_m^T     f(n)     (b'_m)^T 
      H_m^T     b'_m     A_m + f(n)I_m ,

where b_m=m^2     (m-1)^2    …    1^T, b'_m=K_mb_m and H_m=(h_ij) ∈ M_m() is defined by
h_ij=(m+1+j-i)^2. We will write U as: 
 
    UU^T=VV^T+V'V'^T+tt^T+ZZ^T,

where (as in the proof of Lemma <ref>) each column of V will correspond to an off-diagonal pair of entries in the upper-left block of B_n, and each column of V' will correspond to an off-diagonal pair of entries in the lower-right block of B_n. In this case we will need to consider separately the position (m-1, m) in the upper-left block, and  the position (m+2, m+3) in the lower-right block. To deal with those positions we introduce a rank one matrix tt^T. Finally,  each column of Z will correspond to an entry in the upper-middle block of B_n.

As in the even case, we define
V:=∑_i=1^m-2∑_j=i+1^m u_ij u_ij^T,
where 

    u_ij   := (j-i) (e_i + e_j + α_ij e_m+1+k),  
     k:= ⌊m+j+1/2⌋ - i, 
    α_ij   := 1/k(2m-i-j+2),

and V':=K_nV.
There are two positions in the diagonal blocks left to cover, and this is done by tt^T, where:
t := e_m-1+e_m+e_m+2+e_m+3. It is easy to verify that V^Tw(n)=0, V'^Tw(n)=0, t^Tw(n)=0, and

    VV^T+V'V'^T+tt^T= A_m + Ê_1     0    Ŝ
     0     0     0 
    Ŝ    0     A_m+Ê_2 ,

where Ê_1 and Ê_2 are diagonal matrices. Finally, we define 
Z:=∑_i=1^m z_i z_i^T, where

    z_i := (m+1-i) (e_i + e_m+1 + e_2m+2-i).

We have z_i^Tw_n=0, and 

    ZZ^T = F     b_m     FK_m 
     b_m^T    β    (b'_m)^T  
     K_mF     b'_m     K_mFK_m ,

where F is a diagonal matrix.

To complete the proof, we need to show that α:= f(n)- β > 0, and that the matrix C:=H_m-Ŝ-FK_m is non-negative. 
We have 

    β = ∑_i=1^m i^2 = 1/6m(m+1)(2m+1) < f(n) = 2/3m(m+1)(2m+1),

proving α >0. 









The matrix tt^T influences only the lower-left 2 × 2 corner of C. For n=5 we have

    C= 8     11 
     2     8,
 and for n ≥ 7 the lower-left 2 × 2 corner of C is equal to 
    0     6 
     2     0≥ 0.



The (x,y) entry of the matrix FK_m is nonzero only when x+y=m+1, and in this case (FK_m)_xy = y^2.

We will consider VV^T and V'V'^T separately, and to this end we denote by S the top-right corner of VV^T and by S'=K_mS^TK_m the top-right corner of V'V'^T. Contributions of S and S' are dealt with as in Lemma <ref>, with only minor changes in indexing. 


Considering nonzero entries of the matrix u_iju_ij^T at position (x, m+1+y), we get (x, y) = (i, k) or (x, y) = (j, k). This results in three possibilities, as follows.  
First we can have i=x and j= 2x+2y-1-m, and we define the corresponding matrix S_I by 

    (S_I)_xy=(u_x,2x+2y-1-mu_x,2x+2y-1-m^T)_x,m+1+y
 for 1 ≤ x < 2x+2y-1-m ≤ m, x≤ m-2,
and zero for (x, y) outside this  range. Similarly, for  i=x and  j= 2x+2y-m we have S_II with nonzero entries  

    (S_II)_xy=(u_x,2x+2y-mu_x,2x+2y-m^T)_x,m+1+y
 for 1 ≤ x < 2x+2y-m ≤ m, x≤ m-2, and for  j=x and i = ⌊m+x+1/2⌋ - y we define S_III with nonzero entries 

    (S_III)_xy=(u_⌊m+x+1/2⌋-y,xu_⌊m+x+1/2⌋-y,x^T)_x,m+1+y
  for 1 ≤⌊m+x+1/2⌋-y < x ≤ m, ⌊m+x+1/2⌋-y ≤ m-2.











Computation of relevant entries of u_iju_ij^T and reordering of the bounds, gives us the following formulas for the nonzero entries of S_I, S_II and S_III:

    (S_I)_xy   = 1/y(x+2y-m-1)^2(3m-3x-2y+3)  for   m+2-2y ≤ x ≤ m-y, 
    
    (S_II)_xy   = 1/y(x+2y-m)^2(3m-3x-2y+2)  for   m+1-2y ≤ x ≤ m-y, y ≥ 2, 
     
    (S_III)_xy   =1/y(x-⌊m+x+1/2⌋+y)^2(2m-⌊m+x+1/2⌋+y-x+2) 
        for max{m+2-2y,2y-m}≤ x ≤ m, (x, y)(m, 1).


We estimate the matrix S_III by a matrix Ŝ_III, Ŝ_III≥ S_III, with nonzero entries in positions (x,y) satisfying m+2-2y ≤ x ≤ m, (x, y)(m, 1) and defined by: 

    (S_III)_xy≤ (Ŝ_III)_xy = 1/y(x-m+x/2+y)^2(2m-m+x/2+y-x+2).


The matrices H_m, S + S' and FK_m are symmetric with respect to the counter-diagonal, so it is enough to prove the nonnegativity of elements (H_m-FK_m-S-S')_xy that satisfy x+y ≤ m+1. Inside this region we will consider four subregions, that are determined by the location of nonzero entries of matrices S_I, S_II, S_III, S'_I, S'_II, S'_III and FK_m. Since (S'_I)_xy=(S'_II)_xy=0 for x+y≤ m+1, matrices S'_I and S'_II don't need to be considered. The four regions are defined as follows:
 

  * {(x,y);  x ≤ m-2y}, only (S'_III)_xy≠ 0, hence C_xy=(H_m-S'_III)_xy≥ (H_m-Ŝ'_III)_xy.  Introducing new variables u=y-1≥ 0 and v=m-x-2y≥0 gives us:

    C_xy≥24 + 220 u + 258 u^2 + 81 u^3 + 8 v + 104 u v + 66 u^2 v + 12 u v^2/8 (3 + 2 u + v) > 0.


  * {(x,y);  x = m+1-2y, y ≥ 2}, (S_II)_xy≠ 0 and (S'_III)_xy≠ 0, hence C_xy≥ (H_m-S_II-Ŝ'_III)_xy, which we express in terms of u=y-2≥ 0. We get

    C_xy≥305 + 691 u + 435 u^2 + 81 u^3/16 (2 + u) > 0.


  * {(x,y);  m+2-2y ≤ x ≤ m-y}, (S_I)_xy≠ 0, (S_II)_xy≠ 0, (S_III)_xy≠ 0, (S'_III)_xy≠ 0, and C_xy≥ (H_m-S_I-S_II-Ŝ_III-Ŝ'_III)_xy.
 Change of variables u=x+2y-m-2≥ 0 and v=m-x-y≥0 gives us

    C_xy≥p(u,v)/8 (2 + u + v) (3 + u + 2 v) ,

where

    p(u,v)    = -122 - 71 u + 25 u^2 + 30 u^3 + 6 u^4 + 361 v + 568 u v + 241 u^2 v  + 
       +  36 u^3 v + 909 v^2 + 825 u v^2 + 156 u^2 v^2 + 531 v^3 + 258 u v^3 +  81 v^4.

We get C_xy > 0 unless (u, v) ∈{ (0, 0), (1, 0)}. For (u, v) ∈{ (0, 0), (1, 0)} we get (x, y) ∈{(m-2,2), (m-3, 3)}. In those two cases we compute C_m-2,2 = 3 >0 and C_m-3,3 = 11/4 >0.

  * {(x,y);   x=m+1-y, y ≥ 3 }, (FK_m)_xy≠ 0, (S_III)_xy=(S'_III)_xy≠ 0, and C_xy≥ (H_m-FK_m-2 Ŝ_III)_xy, 
which we express in terms of u=y-3≥ 0 to get

    C_xy≥68 + 116 u + 52 u^2 + 7 u^3/4 (3 + u) > 0.






The matrix B_n = A_n + f(n)I is completely positive for every n.




Suppose first that n=2m is even. By Lemma <ref> we write B_n as 

    B_n = UU^T + D     C 
     C^T     D',

and since the matrix D     C 
 C^T     D' satisfies conditions of Lemma <ref>, the proof is complete in this case. 

Suppose now that n=2m+1. By Lemma <ref> we can write B_n as 

    B_n = UU^T + D     0     C 
     0    α    0 
     C^T     0     D'.

Now D     C 
 C^T     D' is completely positive by Lemma <ref>, since its eigenvector corresponding to 0 is just w(2m+1) with the middle entry (which is equal to zero) omitted. In addition, we have proved that α > 0, and this allows us to conclude that

    D     0     C 
     0    α    0 
     C^T     0     D'
 is completely positive.









§  INTEGER COMPLETELY POSITIVE FACTORIZATION 


In this section we want to take into account that A_n has integer entries, and we ask the question, for what values of g(n) does the matrix A_n + g(n)I have an integer completely positive factorization. It turns out, that for small n,  g(n)=f(n) works, but for n≥ 6 we need to choose g(n) >f(n).


Let B_n:=A_n+f(n)I_n. Then:


  * B_n has an integer completely positive factorization for n ≤ 5.

  * B_6 does not have an integer completely positive factorization.

  * B_6+I_6 has an integer completely positive factorization.




To prove the first item and the third item, we find integer completely positive factorizations explicitly, as follows: 

    B_2    =  1     1 
     1     1  =  1 
     1  1     1 ,

 
    B_3    =  4     1     4 
     1     4     1 
     4     1     4  = 
     1 
     1 
     1  1     1     1  + 
    3  1     0 
     0     1 
     1     0  1     0     1 
     0     1     0 ,

 
    B_4    =  10     1     4     9 
     1     10     1     4 
     4     1     10     1 
     9     4     1     10  = 
     1     1     0 
     1     0     3 
     1     3     0 
     1     0     1  1     1     1     1 
     1     0     3     0 
     0     3     0     1  + 
    8  1 
     0 
     0 
     1  1     0     0     1  ,

 
    B_5    =  
    20     1     4     9     16 
     1     20     1     4     9 
     4     1     20     1     4 
     9     4     1     20     1 
     16     9     4     1     20  = 
     1     2     0     0 
     1     0     4     0 
     1     0     0     4 
     1     4     0     0 
     1     0     2     0  1     1     1     1     1 
     2     0     0     4     0 
     0     4     0     0     2 
     0     0     4     0     0  + 
       + 3  1     0     2 
     0     1     0 
     1     0     0 
     0     1     0 
     1     0     2  1     0     1     0     1 
     0     1     0     1     0 
     2     0     0     0     2  .

Furthermore: 
 
    B_6+I_6 =A_6+36 I_6=
      
    36     1     4     9     16     25 
     1     36     1     4     9     16 
     4     1     36     1     4     9 
     9     4     1     36     1     4 
     16     9     4     1     36     1 
     25     16     9     4     1     36  
     =
       =  1 
     1 
     1 
     1 
     1 
     1  1     1     1     1     1     1  + 
    2  1     0     0 
     0     0     2 
     0     4     0 
     4     0     0 
     0     0     2 
     0     1     0  1     0     0     4     0     0 
     0     0     4     0     0     1 
     0     2     0     0     2     0  + 
       + 3  1     0 
     0     1 
     1     0 
     0     1 
     1     0 
     0     1  1     0     1     0     1     0 
     0     1     0     1     0     1  +
    6  1     0     2 
     0     2     0 
     0     0     0 
     0     0     0 
     2     0     0 
     0     1     2  1     0     0     0     2     0 
     0     2     0     0     0     1 
     2     0     0     0     0     2  .


To prove the second item, we assume that  

    B_6=  
    35     1     4     9     16     25 
     1     35     1     4     9     16 
     4     1     35     1     4     9 
     9     4     1     35     1     4 
     16     9     4     1     35     1 
     25     16     9     4     1     35

has an integer completely positive factorization: B_6 = UU^T. Let u_i, i=1,…,r, denote the columns of U. First we note, that all columns of U need to be orthogonal to
w =  5     3     1     -1     -3     -5 ^T by Proposition <ref>.  
Next we look at conditions that are coming from the fact that the superdiagonal of B_6 has all elements equal to 1. In particular, this implies, that if two consecutive entries of any column of U are nonzero, they both have to be equal to 1.  
Moreover, since (B_6)_12 = 1, one 
of the columns in U has to have the first two entries both equal to one, and, without loss of generality,  we may assume that this holds for the first column of U: 
u_1 =  1     1     *     *     *     * ^T . Taking into account all the conditions on u_1 that we have listed so far, u_1 has to be equal to one of the following three vectors: 

    v    =  1     1     1     1     1     1 ^T ,
    
    v'    =  1     1     0     0     1     1 ^T , 
    
    v”   =  1     1     0     3     0     1 ^T .

In the first case we have 

    B_6 - u_1u_1^T = U_1U_1^T =  
    34     0     3     8     15     24 
     0     34     0     3     8     15 
     3     0     34     0     3     8 
     8     3     0     34     0     3 
     15     8     3     0     34     0 
    
    24     15     8     3     0     34
     .


If u_1=v', then some column of U (suppose u_2) has 1 in the second and the third position. This forces the first entry of u_2 to be zero, and the condition w^Tu_2=0 leaves us with only one choice for u_2: u_2 =  0     1     1     1     1     0 ^T , thus 

    B_6 - u_1u_1^T - u_2u_2^T = U_2U_2^T =  
    34     0     3     8     15     24 
     0     33     0     2     7     15 
     3     0     33     0     2     8 
     8     2     0     33     0     3 
     15     7     2     0     33     0 
    
    24     15     8     3     0     34
     .


If u_1=v”, we notice that some column of U (suppose u_2) has 1 in the last two positions, and since the first two entries of u_2 cannot both be equal to one, we deduce that
u_2 =  1     0     3     0     1     1 ^T . Since (B_6)_23=1, we still need a column, say u_3, in U that has 1 in both the second and the third position. Again, we are left with only one option:
u_3 =  0     1     1     1     1     0 ^T , and

    B_6 - u_1u_1^T - u_2u_2^T - u_3u_3^T = U_3U_3^T =  
    33     0     1     6     15     23 
     0     33     0     0     8     15 
     1     0     25     0     0     6 
     6     0     0     25     0     1 
     15     8     0     0     33     0 
    
    23     15     6     1     0     33
     .



In all three cases we are left with a matrix U_i U_i^T that has all the entries on the superdiagonal equal to zero. In particular, 
the columns of U_i have at least every other entry equal to zero. Moreover, one of the columns in U_i has to have the first and the third
entry nonzero, so it is of the form:
v̂=  a     0     b     0     c     d ^T ,
where cd=0. 
Looking at the values of U_iU_i^T in positions (1,3), (3,5) and (2,6) we get ab≤ 3, bc ≤ 3, and  bd ≤ 8. This, together with the condition  w^Tv̂=0, gives us only one option for v̂:

    v̂ =  1     0     1     0     2     0 ^T .

In the third case we also need bc=0, so this choice immediately leads to contradiction. In the remaining two cases U_iU_i^T - v̂v̂^T still has the (1,3)-entry positive, so we need to have another column in U_i, i=1,2,  of the same form. 
However, in both cases U_iU_i^T-2 vv^T has a negative entry in (3,5) position, leading to a contradiction. 


The ad hoc approach, that we used to find completely positive factorizations in the proof above, is not well suited for generalisation to matrices of arbitrary size. Instead, we develop a more general technique, that can be applied to a class of Toeplitz matrices, but that does not produce the optimal q for distance matrices. 


Let n ∈ℕ, i ∈{1,…,n-1}, and 
    E_i := I + J_n^i + J_n^2i + ... + J_n^(n-1)i + (J_n^T)^i + (J_n^T)^2i + ... + (J_n^T)^(n-1)i .

Any matrix of the form ∑_i=1^n-1 a_i E_i, a_i ∈ℕ∪{0}, has an integer completely positive factorization. 



Clearly, it is sufficient to prove the statement for E_i, i=1,…,n-1. So,  let us choose i ∈{1,…,n-1}, and write
n = qi+r where r ∈{0,…,i-1}. We define: 

    U_i^T =  I_i     I_i     ...     I_i     D_r ,

where D_r is an i × r (possibly null) matrix with 1's on the main diagonal and 0's elsewhere. Equality E_i=U_iU_i^T can be checked by a straightforward calculation. 


Before we see how this lemma can be applied to distance matrices, we illustrate its application by an example.


The matrix of the form 
    E_1 + E_2 + … + E_n-1 = (n-1)I + ∑_i=1^n-1τ(i)(J_n^i + (J_n^T)^i)

has an integer completely positive factorization, where 
τ(i) := ∑_d | i d^0, the number of all divisors of i.








The previous example shows a straightforward application of Lemma <ref>, where a_i are all chosen to be one. To determine the values for a_i, i=1,…,n-1, so that the matrices A_n and ∑_i=1^n-1 a_i E_i agree outside the diagonal, we need the Jordan totient function: 

    J_2(k) :=k^2 ∏_p | k 
     p  prime(1-1/p^2) .



The matrix A_n + g_J(n)I has an integer completely positive factorization for g_J(n):= ∑_k=1^n-1 J_2(k).



The claim will follow from Lemma <ref>, after we show that:
 
    A_n + g_J(n)I_n=∑_i=1^n-1J_2(i)E_i.

Clearly, the equality holds for the elements on the main diagonal. The equality for off-diagonal elements follows from the following well known formula, that can be found for example in <cit.>: 

    ∑_d  divisor of  k J_2(d) = k^2.






We know that, g_J(n) is not optimal, since, for example, g_J(6)=48>36, and 36 is the optimal choice for n=6 by Proposition <ref>. 
On the other hand g_J(n) <g_D(n), where g_D(n):=1/6n(n-1)(2n-1) is the smallest value that makes A_n+g_D(n)I_n diagonally dominant.  Asymptotically we have 

    g_J(n)∼n^3/3ζ(3),
 
where ζ is Riemann zeta function, see for example <cit.>.
This gives us: 
    g_J(n) ∼2/ζ(3)f(n) ≈ 1.66381f(n),

which is again an improvement over g_D(n) ∼ 2 f(n). 

amsplain












