2
$\begingroup$

Which method should I choose to solve a large (~20 000 variables) dense symmetric positive-definite, possibly ill-conditioned, system of linear equations? The system will be solved for two vectors. I'm interested both in performance and numerical accuracy, so I'm looking for possible solutions to choose the best compromise. I want to use it as a part of a special purpose SDP solver under OCTAVE/MATLAB, so most probably it will be implemented as a mex file, with a single thread.

EDIT: Thank you for your responses. I mean an optimization solver, with primal-dual interior point method. Dual solution is always feasible in my case, and dense (I mean the dual is not block-diagonal).

EDIT2: The reason why I'm writing an SDP solver is that I want to try to get a more efficient tool for solving large NPA problems (http://arxiv.org/abs/0803.4290). This case allows for some simplifications, the interior is always non-empty, I can find a good initial dual feasible solution (and thus keep this feasibility), and the structure allows for easier creation of the equation for dy. The problem is that this equation is large for these problems.

I tried to use SeDuMi and SDPT3, with the latter being able to run much larger programs. For large I got with SDPT3 the error saying "Schur complement matrix not positive definite" meaning that Cholesky factorization failed.

EDIT3: I observe that with this warm start it takes about 3-4 iterations less to solve the program. The problem is the failure of Cholesky for the large Schur complement matrix. In this case it is fully dense, but near the last iteration it becomes very bad conditioned. This is why I'm looking for a method which will be able to solve it.

EDIT4: Thank you for your response. I'll try pivoting with Cholesky. Do you recomend any package which works under OCTAVE for this task?

Improve this question
$\endgroup$
8
  • $\begingroup$ Knowing more about the special purpose SDP solver might help in answering this question. Are you implementing a primal-dual interior point method? Or, a first order method that requires projections onto the set of solutions to $A(X)=b$? Or something else? $\endgroup$
    – Brian Borchers
    Commented Apr 18, 2015 at 18:52
  • $\begingroup$ I think the OP meant SPD (for the matrix), not SDP (for the problem where it comes from). $\endgroup$
    – Wolfgang Bangerth
    Commented Apr 19, 2015 at 18:32
  • 1
    $\begingroup$ Perhaps the OP will edit his question to provide more information... $\endgroup$
    – Brian Borchers
    Commented Apr 19, 2015 at 18:37
  • $\begingroup$ In response to "EDIT2" Note that having a dual feasible solution won't be particularly helpful if that solution is not well centered. In general, the primal-dual interior point method doesn't "warm start" well. Having problem structure that makes it easy to construct the Schur complement matrix could potentially be helpful. $\endgroup$
    – Brian Borchers
    Commented Apr 23, 2015 at 14:47
  • 1
    $\begingroup$ A final trick that can help immensely is to determine the pivot that results in failure of the Cholesky factorization and simply force that element of $\Delta y$ to be 0 (e.g. by making the corresponding diagonal element very large.) $\endgroup$
    – Brian Borchers
    Commented Apr 23, 2015 at 23:30

2 Answers 2

Reset to default
6
$\begingroup$

In double precision, a system of 20,000 equations in 20,000 variables will require 3.2 gigabytes of RAM to store the matrix. This isn't "big" by contemporary standards, and is reasonably easy to solve in practice.

To make this run fast, you'll want to use of a well tuned linear algebra library. The standard for this kind of work would be to use LAPACK and BLAS routines, with a well tuned, multithreaded, implementation of the BLAS. You'd use the Cholesky factorization routine in LAPACK to compute a Cholesky factorization and then use other LAPACK routines to do the forward and backward solves.

In the MATLAB environment, your best bet is to piggy-back off of MATLAB. The Mathworks supplies efficient BLAS/LAPACK routines with MATLAB and handles any licensing issues. You can call these routines from within your MEX routine.

In Octave you should hope that the installation of Octave came with a good multithreaded BLAS. Although they aren't quite as good as Intel's MKL or AMD's ACML libraries (which are not open source), there are open source alternatives including Openblas and ATLAS.

If you're implementing a primal-dual interior point method for SDP, you should be aware that some special tricks are needed to make Cholesky factorization work well. You should also be aware that in many cases, the system of equations actually isn't dense, and it can be worthwhile to use a sparse Cholesky factorization routine.

Improve this answer
$\endgroup$
2
  • 1
    $\begingroup$ Since you've said that you're implementing a primal-dual interior point method for SDP, a couple of comments. The large symmetric positive definite system of equations tha arises in the primal-dual interior point method is sometimes sparse, and it can be very worth while to use a sparse Cholesky factorization routine in those situations. Otherwise, BLAS/LAPACK is where it's at. Iterative methods have not worked well in this situation, because of the extreme ill-conditioning of the Schur complement matrix. $\endgroup$
    – Brian Borchers
    Commented Apr 23, 2015 at 0:35
  • $\begingroup$ A more obvious question is why you aren't making use of one of the existing open source SDP solvers? $\endgroup$
    – Brian Borchers
    Commented Apr 23, 2015 at 0:38
1
$\begingroup$

Provided the matrix is symmetric and positive-definite, the conjugate gradient method is the provably fastest iterative method, requiring $O(\sqrt\kappa)$ iterations to converge. Depending on how badly-conditioned your matrix is, that may be more favorable than using a direct method like the Cholesky decomposition to solve your system. It is likely that you'll need some kind of preconditioner, the choice of which will depend on the nature of your problem.

A 20,000 x 20,000 dense matrix is getting uncomfortably close to the limit of what can be addressed on a 32-bit machine. Storing that matrix alone will require 3.2 GB of memory. All that you need for many iterative methods is to be able to compute the matrix-vector product $A\cdot x$; for some problems, this can be done without explicitly creating the matrix $A$ entry-by-entry. If you were able to do that for your problem, I think it would speed things up dramatically.

Failing that, I think you'll need to resort to multiple threads and an optimized BLAS/LAPACK implementation like ATLAS if you want good performance. You may also want to look into the library Elemental, which does distributed dense linear algebra and has a C interface.

Improve this answer
$\endgroup$
3
  • 2
    $\begingroup$ It's not true that CG is the provably best method for SPD matrices. This lives off the assumption that you can do everything in around $\frac 13 N$ iterations without a preconditioner for a dense matrix. Otherwise a simple Cholesky decomposition will be faster. If you want to employ a preconditioner, assuming the preconditioner takes only as much time as a matrix-vector product (that's about the cheapest it can get) then you need to use fewer than $\frac 16 N$ iterations, and even less if your preconditioner is more expensive. $\endgroup$
    – Wolfgang Bangerth
    Commented Apr 19, 2015 at 18:36
  • $\begingroup$ I meant provably fastest Krylov subspace method, not that it's faster than Cholesky. $\endgroup$
    – Daniel Shapero
    Commented Apr 20, 2015 at 1:19
  • $\begingroup$ Even that I would doubt. How would you prove that it's the fastest method? $\endgroup$
    – Wolfgang Bangerth
    Commented Apr 20, 2015 at 1:52

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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