The parametric primal-dual simplex algorithm

In a recent Subject To podcast featuring Bob Vanderbei, then Bob mentions a primal-dual simplex algorithm. Bob also discusses the algorithm in his Linear programming book and has been teaching it. It is clear from the podcast that he really likes this algorithm.  

Now I have recently also gotten interested in that algorithm for the following reasons:

• It is a single phase algorithm that is trivial to start. 
• At least in theory it has no problems with degeneracy given an easy to satisfy condition holds.
• It has a convergence measure that monotonically decreases to zero. 
• Given a nondegeneracy assumption is satisfied then the primal basic variables and the dual nonbasic variables are strictly positive i.e. are in the interior. 
• Moreover, the primal and dual infeasibility measures are reduced at the same rate.
• Therefore, in some sense it seems closely related to the homogenous model that is employed in interior-point methods.

In version 1 of Bob's LP book, Bob attributes the algorithm to G. B. Dantzig and references his famous 1963 book.  Now Bob also mentions that that the algorithm is a special case of Lemke's algorithm for linear complementarity problems published in 1965. This fact has been proved by Irv Lustig as mentioned in his ph.d. thesis from 1987.

I have been searching for computational results for the primal-dual simplex algorithm but have not located much beyond what is presented by Vanderbei. One exception is the previously mentioned ph.d. thesis of Lustig that contains some of the very limited computational experiences with the primal-dual simplex algorithm. Note how few and small problems Lustig reports results back in 1987. 

Nimrod Megiddo has an interesting paper showing that many seemingly different algorithms can be seen as a special case of the primal-dual simplex algorithm. The basic idea is that by choosing the initial perturbation employed by the primal-dial simplex algorithm appropriately the behavior of the other algorithms can be recovered.  I doubt that for instance the dual simplex algorithm using steepest-edge pricing is  a special case though.

I shall not claim that the primal-dual simplex algorithm is better in practice than the well known primal and dual simplex algorithms. However, from a theoretical point of view I find the primal-dual simplex algorithm much more pleasing.   

I may add more material to this blog post later.

Matrix ordering and graph partitioning

Since I spend a lot of time on solving sparse linear equation systems then I am also a user of sparse matrix reordering methods. My claim to fame is that I have implemented approximate minimum degree myself and it is used in MOSEK.

Below I summarize some interesting link to graph partitioning software:

• Metis (the most famous).
• Mondriaan (do not know much about this).
• Spooles (Not that active).
• Multisection a generalization of bisection.

My Git links and tips

At my company we have started using the source control system Git instead of Perforce. Here I have collected some Git related links.

• Understanding Git conceptually

Installing TexLive 2016 on Windows Server 2016

My first attempt at intsalling TexLive 2016 on a Windows 2016 failed. However, the following worked

• Download install-tl-windows.exe from http://tug.org/texlive/acquire-netinstall.html
• set PATH=c:\windows\system32;%PATH%
• Run the script install-tl-windows.bat

then reason is that cmd.exe has to be on the PATH.

Error handling in BLAS libraries

It is very common to use a BLAS library to perform linear algebra operations such as dense matrix times dense matrix multiplication which can be performed using the dgemm function. The advantage of BLAS is

• it is well a defined standard.
• and hardware vendors such as Intel supplies a tuned version.

Now at MOSEK my employeer we use the Intel MKL library that includes a BLAS implementation. It really helps us deliver good floating point performance. Indeed we use a sequential version of Intel MKL but call it from potentially many threads using Clik plus. This works well due to the well designed BLAS interface. However, there is one rotten apple in the basket and that is error handling.

Here I will summarize why the error handling in the BLAS standard is awful from my perspective.

First of all why can errors occur when you do the dgemm operation if we assume the dimensions of the matrices are correct and ignoring issues with NANs and the like. Well, in order to obtain good performance the dgemm function may allocate additional memory to store smallish matrices that fit into the cache. I.e. the library use a blocked version to improve the performance.

Oh wait that means it can run of memory and then what? The BLAS standard error handling is to print a message to stderr or something along that line.

Recall that dgemm is embedded deep inside MOSEK which might be embedded deep inside a third party program. This implies an error message printed to stderr does not make sense to the user. Also the user would NOT like us to terminate the application with a fatal error. Rather we want to know that an out of space situation happened and terminate gracefully. Or doing something to lower the space requirement. E.g. use a fewer threads.

What is the solution to this problem? The only solution offered is to replace a function named xerbla that gets called when an error happens. The idea is that the function can set a global flag indicating an error happened. This might be a reasonable solution if the program is single threaded. Now instead assume you use a single threaded dgemm (from say Intel MKL) but call it from many threads. Then first of all you have to introduce a lock (a mutex) around the global error flag leading to performance issues. Next it is hard to figure out which of all the dgemm calls that failed. Hence, you have to fail them all. What pain.

Why is the error handling so primitive in BLAS libraries. I think the reasons are:

• BLAS is an old Fortran based standard. 
• For many years BLAS routine would not allocate storage. Hence, dgemm would never fail unless the dimensions where wrong.
• BLAS is proposed by academics which does not care so much about error handling. I mean if you run out of memory you just buy a bigger supercomputer and rerun your computations.

If the BLAS had been invented today it would have been designed in C most likely and then all functions would have returned an error code. I know dealing with error codes is a pain too but that would have made error reporting much easier for those who wanted to do it properly.

Here are some links with background information:

• 2009 Intel MKL foum post (The 2 last post is the most interesting)
• 2012 Intel MKL forum post (By me on my problems with error handling)
• 2015 Intel MKL forum post (Some understand on why things as they are)
• BLIS discussion (Newer BLAS library does not focus on error handling either)

Nonsymmetric semidefinite optimization problems.

In  semidefinite optimization we optimize over a matrix variable that must be symmetric and positive semidefinite.

Assume we want to relax the assumption about symmetry. Is that an important generalization? The answer is no for the following reason. Since 

   (X+X')/2 is PSD

implies

   X 

is PSD. Observe

   X = (X+X')/2+(X-X')/2

and

  y'( (X-X')/2) y >= 0.

implying X is PSD.

Note  (X-X')'=-(X-X') implying X-X' is skew symmetric.

Hence, any nonsymmetric semidefinite optimization problem can easily be posed as a standard symmetric semidefinite optimization problem.

Multi threaded programming

I found the talk: Plain Threads are the GOTO of todays computing by Hartmut Kaiser very interesting because I have been working on improving the multithreaded code in MOSEK recently. And is also thinking how MOSEK should deal with all the cores in the CPUs in the future. I agree with Hartmut something else than plain threads is needed.

Here are some potential replacements:

• Berkeley Unified Parallel C
• CheckedThreads  Source with a liberal license. Cilk like. Both C and C++.
• Cilk Plus Supported by Intel C, gcc and Clang.
• HPX C++ only.
• Occa. A unified approach to multithreaded languages. See also.
• OpenMP Supported by many compilers.
• Pfunc Cilk like features. Liberal license. Can be called from C but is C++.
• Starpu
• Threaded building blocks. C++ only.
• Wool  A pure C implementation that provides something very close to Cilk. 
• XKAAPI

Previously I have used OpenMP but I really dislike that. In my opinion it is somewhat ugly and you feel the Fortran origins. Recently I have played with Cilk which is very simple.

Checkedthreads seems like a good option if a simple C only tool can do the job. I have plan to try that at MOSEK.

Pfunc seems very powerful and there is ph.d. thesis about its design. The project might be semi-dead though.

Wool also seems very interesting. It is plain C and the author claims the spawn overhead is very low. There is also an older comparison with other task libraries.

Btw I biased towards tools that has no C++ bindings because currently MOSEK is plain C code. Adding a dependency on a C++ runtime library adds headaches.

Some common terminology when working on parallism is

• Fork-join parallel model
• Pthreads programming and tutorial