Is there a faster way to check whether $A\in\mathbb{R}^{n\times n}$ is positive definite on $b^{\bot}:=\{x\in \mathbb{R}^{n}: x\cdot b=0\}$ than

function foo(A,b)
    N = null(b);
    e = eig(N'*A*N);
    return all(e>0);

Instead of e = eig(N'*A*N);, you can use [R,p]=chol(N'*A*N);, and test for p==0. Matlab function chol returns the cholesky factor R and a p which is zero if N'*A*N is positive definite. p is a positive integer if the matrix is not positive definite.

Usually chol is much faster than eig. For example, with dimension n=2000, I found that chol(N'*A*N) is about 10 times faster than eig(N'*A*N) (tested with Octave instead of Matlab).

The matrix product N'*A*N is relatively expensive to compute. It might be faster to use a sparse basis N for $b^{\bot}$: Matlab/Octave code:

function r = foo2(A,b)
    [dummy,p] = sort(abs(b),"descend");
    A = A(p,p);
    m = size(A,1);
    b = b(p);
    N = [b(2:m)/-b(1);speye(m-1,m-1)];   % now b*N = 0
    [L,p] = chol(N'*A*N);
    r = (p==0);
  • $\begingroup$ That construction of N is brilliant. The sparsity is great of course, but I wasn't even aware about that explicit construction. $\endgroup$ – Bananach Dec 13 '18 at 9:31
  • $\begingroup$ (Now that I think of it, it's kind of a super simple version of Gram-Schmidt starting from $(b,e_2,\dots,e_m)$ without the need for an orthogonality nor normality. Still, awesome) $\endgroup$ – Bananach Dec 13 '18 at 9:34
  • $\begingroup$ Note that N = null(b); produces an orthogonal N (in contrast to the one in foo2), which should be preferred if your original matrix A is ill conditioned. I haven't looked into the details, but I guess that the matrix N = [b(2:m)/-b(1);speye(m-1,m-1)]; has a condition number $\kappa(N) \le \sqrt{m}$. $\endgroup$ – wim Dec 13 '18 at 13:49
  • $\begingroup$ For $b \neq 0$ you can easily check the condition number of the sparse N with cond_N = norm(N*b(2:m)'/norm(b(2:m))), and decide to use N = null(b); if cond_N is too large. $\endgroup$ – wim Dec 13 '18 at 14:21
  • $\begingroup$ (Note that if the division in cond_N = norm(N*b(2:m)'/norm(b(2:m))) fails due to norm(b(2:m)) == 0, then cond_N = 1.) $\endgroup$ – wim Dec 13 '18 at 14:47

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.