I know that the KrylovKit.eigsolve Julia function can find a vector @@x@@ to maximize or minimize @@x^\dagger Ax@@ for a fixed matrix @@M@@. We can think of this in index notation as @@M_{\alpha\beta}x^\alpha\overline{x}^\beta@@, where the bar denotes complex conjugation.

What I am wondering is: does some common package in Julia have a built-in generalization of this method that optimizes over multiple input vectors? For instance, is there a built-in function that finds vectors @@x@@ and @@y@@ to maximize or minimize @@M_{\alpha\beta\gamma\delta}x^\alpha\overline{x}^\beta y^\gamma\overline{y}^\delta@@?

In case you're wondering why I'm asking, basically I am wondering whether it is possible to do a single optimization step over two non-adjacent sites in a way that keeps them unliked, because if you just multiply the two tensors together and then optimize that, then it will introduce a link between those sites. Basically I want to just do this step once in a while to see whether it is possible to fix some problems with my trial Schwinger model (as I have already discussed in some other questions on this website).