How to solve the linear system in your TT-IPM?

[optztgit]

Very interesting work :)

But still confused about the linear system solver. Is it a direct AMEn procedure [1, algorithm 4] with iteratively updating each tensor-train block of W?
(It seems that the linear system in TT-IMP is not SPD as required in [1].) If not, how to solve it with the TTM structure and change the TT rank during the iterations?

Looking forward to your reply!

Ref: [1]Dolgov, S., & Savostyanov, D. (2014). Alternating minimal energy methods for linear systems in higher dimensions. SIAM Journal on Scientific Computing, 36(5), A2248-A2271.

[FreditorK]

Thanks for the comment!

The solver here is not a direct implementation of Algorithm 4 in [1]. Instead, we use a block-AMEn variant with minor modifications (following [2]) that is better suited for our setting (because of degenerate projectors). A Block-AMEn projects each TT block separately, then reassembles the projected blocks into a local block matrix, which is solved by a local solver to compute the solution update for the block TT tensor.

While the original AMEn theory in [1] requires SPD matrices, in practice AMEn usually outperforms TT-GMRES on nonsymmetric systems, such as the Chemical Master Equation example showcased in [1]. The caveat is that one loses the formal convergence guarantees. To obtain convergence guarantees, one would have to opt for a TT-GMRES.

By shifting the block index from core to core, you achieve a similar TT rank adaptation as in DMRG algorithms. Additionally, we also have the residual enrichment, as in AMEn, that contributes to rank increases before SVD truncates the ranks.

I hope this helps - let me know if it’s clear.

References
[1] S. Dolgov & D. Savostyanov (2014), Alternating minimal energy methods for linear systems in higher dimensions, SIAM J. Sci. Comput. 36(5), A2248–A2271.
[2] P. Benner, S. Dolgov, A. Onwunta, M. Stoll (2020), Low-rank solution of an optimal control problem constrained by random Navier-Stokes equations, Int. J. Numer. Meth. Fluids 92(11), 1653–1678, https://doi.org/10.1002/fld.4843

[optztgit]

Thanks for your very detailed response! I think I have got the key idea of the solver and rank adaption via block TTM index shifting.

I would be grateful if you could help me with two other minor issues:

i) In reference [2], the authors ensure the nonsingularity of the block-reduced system (38) via Theorem 1, which seems to be proved based on some physical background of the matrix parameters (K, K*, M1, M2, and M3). How can the nonsingularity of the k-th small linear system of TT-IPM be ensured?

ii) For the three types of problems mentioned in Section 4 of your paper (max cut, max stable set, and clustering), how do you allocate the indices n_1, ..., n_{2d} for a TTM variable? As discussed at the beginning of [2, Section 4], the allocation is critical and can be related to some physical factors (though the order of the indices is not discussed in [2]).

[FreditorK]

Hi,

Of course, it’s great to see your careful reading of our work; I think I can clarify the two points you mentioned:

i) For the operators $L_Z$ and $L_X$, as they are positive definite, their projections remain positive definite and therefore full rank by the Poincaré separation theorem. So the Cholesky decomposition is never an issue, even in Block-ALS. The other blocks would normally pose an issue, but mostly in the Block-ALS case, i.e. without residual enrichment. Reference [2], to my knowledge, provides the nondegeneracy proofs, as the problem offers this opportunity, and because it is generally an elegant supporting statement. In practice, they also use a residual enrichment, though. For our other blocks, it again comes down to the residual enrichment, which is zero only at the global optimum and otherwise enriches the column space of the projectors, i.e. $\text{span}(U)$ in $R_U$, making degenerate projections, and therefore, convergence to local minima, very unlikely and, in practice, a non-relevant problem.

Additionally, in the code, we include a very small regularisation term on the projected Schur-complement to further ensure nonsingularity. It is a valid point of criticism that this regularisation is not mentioned in the paper, although it is used in practical implementations.

ii) The multiplier allocation is implied by the respective $K_y$ and $K_t$. In maxcut, for example, we only need $n$ multipliers to allocate to the multiplier matrix of size $n \times n$. You could allocate them to the first row of the matrix, or spread them somewhat randomly, or simply position them on the diagonal, which then has implications on the rank of the constraint operators and the $K$'s. One of the best examples is $\mathcal{A}_{part 1}$ of Maximum Stable Set, which is bounded by the TT rank of the graph of the adjacency matrix. If we chose to just take the upper triangular multipliers, as we only need those, the TT-rank is bounded by rank_TT(U A_G) $\leq$ 2 rank_TT(A_G), i.e. it may double the rank, leading to higher rank solutions of the KKT system and also costing more when contracting over. For example, Maxcut has its multipliers on the diagonal, i.e. all indices with i=j in the TT matrix as $K$ is chosen to be vec_TT $(J-I)$.