QR and LQ decompositions

QR decomposition is a best known decomposition from a whole family of orthogonal factorizations, which includes QR, LQ, RQ and QL decompositions. Here Q denotes orthogonal matrix, R and L denote upper and lower triangular matrices. These decompositions may be used to solve full rank least squares problem, as preliminary step in the construction of the SVD decomposition, or for other tasks.

ALGLIB package has C/C++ and C# implementations of two most popular decompositions: QR and LQ. QL and RQ decompositions are less popular and they are not implemented in the ALGLIB.

Aglorithm

QR decomposition functions

Storage format

The most traditional approach to QR decomposition is to represent orthogonal Q as a sequence of Householder reflections. Coefficients of these reflections can be stored in-place, replacing lower triangular entries of input matrix (ones unoccupied by upper triangular R). Well, almost inplace - we need N additional elements denoted Tau. Thus, typical QR decomposition subroutine on input takes square/rectangular matrix A and on output it:

• replaces diagonal and upper triangle of A by triangular factor R
• replaces subdiagonal entries of A by coefficients of Householder reflections which are used to generate/apply orthogonal Q
• returns additional array Tau which stores "addititional" Householder coefficients

ALGLIB functions

QR decomposition functionality is located in ortfac (orthogonal factorizations) subpackage. Following functions can be used with real matrices:

• rmatrixqr - to perform real QR decomposition
• rmatrixqrunpackr - to extract triangular factor R of QR decomposition
• rmatrixqrunpackq - to extract triangular factor Q of QR decomposition
• rmatrixlq - to perform real LQ decomposition
• rmatrixlqunpackl - to extract triangular factor L of LQ decomposition
• rmatrixlqunpackq - to extract triangular factor Q of LQ decomposition

Exactly same set of functions is provided for complex matrices:

• cmatrixqr - to perform complex QR decomposition
• cmatrixqrunpackr - to extract triangular factor R of QR decomposition
• cmatrixqrunpackq - to extract triangular factor Q of QR decomposition
• cmatrixlq - to perform complex LQ decomposition
• cmatrixlqunpackl - to extract triangular factor L of LQ decomposition
• cmatrixlqunpackq - to extract triangular factor Q of LQ decomposition

Benchmarks

Comparison with other free C# libraries

In this section we compare performance of the following free C# libraries:

• ALGLIB 3.12.0 for C# (Free edition) - pure C# version of our product
• Math.NET Numerics 3.20.0 (C# backend) - well known free numerical library, compiled with C# implementation of linear algebra backend
• Accord 3.8.0 - another free numerical library written completely in C#

Although both ALGLIB for C# and Math.NET have high-performance native linear algebra kernels, for the purposes of this benchmark we evaluate performance of pure C# implementations. First, it is interesting to know how much performance NET can bring us. Second, in some cases you have to use 100% managed implementation.

Our first comparison involves single-threaded QR decomposition of 1024x1024 general real matrix. This test was performed on 2.3GHz x64 Intel CPU, running Windows operating system. ALGLIB was compiled with Microsoft C# compiler, for other products we used precompiled NET 4.5 assemblies.

You may see that ALGLIB shows best results close to performance limits of C# language. Accord.NET is ~14x slower, whilst Math.NET is "just" ~3x slower than ALGLIB.

Comparing C#, C++ and SIMD implementations

Another interesting topic to consider is relative performance of C#, generic C/C++ and SIMD-capable implementations. ALGLIB provides several different implementations of linear algebra functionality, all with 100% identical APIs:

• pure C# implementation - 100% .NET, used by ALGLIB for C#
• generic C/C++ one - which can be used on platforms without high-performance SIMD-capable kernels (say, embedded systems)
• HPC core - highly optimized implementation which can be used from C++ and C#, with assembly kernels for performance-intensive parts (including vendor implementations like Intel MKL)

You may guess that these implementations are sorted by performance, with latter being fastest one :) But how fast is it? Our comparison involves single-threaded QR decomposition of 2048x2048 matrix. This test was performed on 2.3GHz x64 Intel CPU, running Windows operating system. Generic C/C++ version is used as baseline one.

You may see that on x64 platform performance of pure C# code is roughly two times lower than that of generic C/C++ code. And, in turn, generic C/C++ code is many times slower than SIMD-capable code utilizing Intel MKL.

This article is licensed for personal use only.