An indefinite LOBPCG type of algorithm for detecting a definite Hermitian matrix pair: MATLAB codes
Description
For a given pair of Hermitian matrices $(A,B)$ of order $n$ with indefinite $B$, this algorithm declares whether this matrix pair is definite or not. A pair of Hermitian matrices is called definite if there exist a real linear combination of these matrices that is a positive definite matrix. Otherwise, it is called indefinite pair. Our algorithm is suitable for medium-size or large and sparse matrix pairs, especially for banded matrix pairs. The proposed algorithm iteratively project a Hermitian matrix pair $(A,B)$ to small dimensional subspaces of $\mathbb{C}^n$ and detect the (in)definiteness of the resulting projected pair. Especially, our algorithm make use of the fact that if the projected pair $(U^HAU, U^H BU)$ is indefinite, then the original pair is also indefinite. The algorithm also include some other conditions to confirm indefiniteness of the given matrix pair or to confirm that the given matrix pair is close to some indefinite matrix pair. On the other hand, the decision that $(A,B)$ is definite is given by finding $\nu\in\R$ such that $A-\nu B$ is a positive definite or negative definite matrix by successfully completed Cholesky factorization. In the algorithm candidates for such $\nu$ are formed from the Ritz values, i.e., from the eigenvalues of the projected pair. The data set contains the MATLAB codes for the algorithm and for the experiments from the paper, accepted for publication in LAA.
Files
Steps to reproduce
-folder pandur_alg has MATLAB files for Algorithm 3.1 - for an initial matrix Xini use a) [Xini,k_pos,k_neg]=InitialGuess(diag(B),k_pos,k_neg) if B=J=diag(+-1) (diagonal matrix of signs) b) [Xini,XiniJ,k_pos,k_neg]=InitialGuessB(B,k_pos,k_neg) if B is not too large (sym_ind_dec(B) is slow) c) [Xini,k_pos,k_neg]=InitialGuess_LDL(B,k_pos,k_neg) if B is large and real symmetric d) some other matrix depending on the structure of B (see test_QEP) - an initial Xini can have more than wanted k_+ and k_- columns - for m=2 in Algorithm 3.1 use [flagP,mu,info]=detecting_definiteness_2dU_dekomp(A,B,Xini,opt,k_pos,k_neg,maxiter,hopr,which_ind_dec,Pcond,tol) - for m=3 in Algorithm 3.1 with B-orthonormalization use [flagP,mu,info]=detecting_definiteness_3dU_dekomp(A,B,Xini,opt,k_pos,k_neg,maxiter,hopr,which_ind_dec,Pcond,tol) - for m=3 in Algorithm 3.1 with B-orthonormalization and the direction of nonpositive curvature use [flagP,mu,info]=detecting_definiteness_3dU_dekomp_NPcurv(A,B,Xini,opt,k_pos,k_neg,maxiter,hopr,which_ind_dec,Pcond,tol) - for m=3 in Algorithm 3.1 with the standard orthonormalization use [flagP,mu,info]=detecting_definiteness_3dU_orth4(A,B,Xini,opt,k_pos,k_neg,maxiter,hopr,which_ind_dec,Pcond,tol) - OUTPUT: flagP = 0 if the pair (A,B) is definite = 1 if the pair (A,B) is not definite = 2 if the algorithm reaches maximal allowed number of iterations = 3 if the pair (A,B) is close to an indefinite one mu = definitizing shift, if pair is definite (a real number for which $A-mu*B$ is a positive definite or a negative definite matrix); [ ] else info = structure information of the algorithm Examples from the paper: - Example 5.1 is in test_clement_indefinite.m - Example 5.2 is in the folder test_QEP - Example 5.3 is in the folder test_saddle_point_form - Example 5.4 is in the folder test_large_banded (bw50.mat includes A and B from the example)