CreateQFromSparseQR.cs

// See https://github.com/wo80/CSparse.NET/wiki/Factorization-members

namespace ConsoleApp
{
    using CSparse;
    using CSparse.Double;
    using CSparse.Double.Factorization;
    using System;
    using System.Diagnostics;
    using System.Globalization;
    using System.Threading;

    static class Test
    {
        public static void Run()
        {
            Thread.CurrentThread.CurrentCulture = CultureInfo.InvariantCulture;

            var A = (SparseMatrix)SparseMatrix.OfArray(new double[,]
            {
                { 7, 3, 1, 0 },
                { 0, 2, 0, 2 },
                { 0, 0, 1, 0 },
                { 1, 0, 0, 2 }
                //{ 0, 1, 1, 1 }
            });

            Run(A, true);
        }

        static void Run(SparseMatrix A, bool print = false)
        {
            // For underdetermined systems, SparseQR factorizes the transposed matrix
            // and the following code won't work.

            if (A.RowCount < A.ColumnCount) throw new Exception("Underdetermined system not supported.");

            var timer = Stopwatch.StartNew();

            PrintMatrix("A:", A);

            var qr = SparseQR.Create(A, ColumnOrdering.MinimumDegreeAtA);

            var t1 = timer.Elapsed;

            qr.GetFactors(out SparseMatrix H, out SparseMatrix R, out double[] beta, out int[] p, out int[] q);

            timer.Restart();

            var Q = CreateQ(H, beta);

            var t2 = timer.Elapsed;

            // Uncomment for benchmark.
            //Console.Write("N = {0,5}: A {1,4:0.0}% , H {2,4:0.0}% , R {3,4:0.0}% , Q {4,4:0.0}%, ", A.RowCount, D(A), D(H), D(R), D(Q));
            //Console.Write("fac time: {0,6:0}ms, Q time: {1,6:0}ms", t1.TotalMilliseconds, t2.TotalMilliseconds);
            //Console.WriteLine();

            if (print)
            {
                var result = Q.Multiply(R.PermuteColumns(Permutation.Invert(p)));

                result.DropZeros(1e-12);
                result.PermuteRows(Permutation.Invert(q));

                PrintMatrix("Q*R:", (SparseMatrix)result);
            }
        }

        static double D(SparseMatrix A)
        {
            return A.NonZerosCount * 100.0 / (A.RowCount * A.ColumnCount);
        }

        static SparseMatrix CreateQ(SparseMatrix H, double[] beta)
        {
            int rows = H.RowCount;
            int cols = H.ColumnCount;

            var hp = H.ColumnPointers;

            int begin = hp[0], end, nnz = 0;

            // Get max nnz of Householder columns.
            for (int i = 0; i < cols; i++)
            {
                end = hp[i + 1];
                nnz = Math.Max(end - begin, nnz);
                begin = end;
            }

            // Dyadic product will have at most nnz^2 entries.
            nnz = nnz * nnz;

            if (!CheckStorageRequirements(rows, cols, nnz))
            {
                Console.WriteLine("Exit.");
            }

            // DyadicProduct assumes row indices to be sorted.
            Helper.SortIndices(H);

            // Allocate all storage outside the loop.
            var E = SparseMatrix.CreateIdentity(rows);
            var D = new SparseMatrix(rows, rows, nnz);
            var P = new SparseMatrix(rows, rows, nnz);
            var Q = (SparseMatrix)E.Clone();

            var Qtemp = new SparseMatrix(rows, rows, nnz);

            // Disable auto-trim for matrix addition and multiplication.
            SparseMatrix.AutoTrimStorage = false;

            for (int i = 0; i < cols; i++)
            {
                DyadicProduct(i, H, D);

                nnz = E.NonZerosCount + D.NonZerosCount;

                if (nnz > Qtemp.Values.Length)
                {
                    Array.Resize(ref Qtemp.RowIndices, 2 * nnz);
                    Array.Resize(ref Qtemp.Values, 2 * nnz);
                }

                E.Add(1.0, -beta[i], D, Qtemp);

                Q.Multiply(Qtemp, P);

                P.CopyTo(Q);
            }

            return Q;

        }

        /// <summary>
        /// 
        /// </summary>
        /// <param name="i">The current column index.</param>
        /// <param name="H">The Householder matrix.</param>
        /// <param name="Q">The Q matrix of the QR factorization.</param>
        /// <remarks>
        /// IMPORTANT:
        /// - the row indices of the Householder matrix <paramref name="H"/> have to be sorted.
        /// - the output matrix <paramref name="Q"/> has to provide enough space to store the non-zero entries.
        /// </remarks>
        static void DyadicProduct(int i, SparseMatrix H, SparseMatrix Q)
        {
            const double EPS = 1e-12;

            // The Householder matrix
            var hp = H.ColumnPointers;
            var hi = H.RowIndices;
            var hx = H.Values;

            // The Q matrix for current column of H.
            var qp = Q.ColumnPointers;
            var qi = Q.RowIndices;
            var qx = Q.Values;

            int columns = Q.ColumnCount;

            // Current non-zero count.
            int nz = 0;

            // Reset the column pointers.
            for (int j = 0; j < columns; j++)
            {
                qp[j] = 0;
            }

            int start = hp[i];
            int end = hp[i + 1];

            // Let v be the i-th column of the Householder matrix (lower triangle).
            // Set the bottom-right submatrix of Q to v*v'.

            // Traverse v'.
            for (int j = start; j < end; j++)
            {
                // Current value of v'.
                double a = hx[j];

                if (Math.Abs(a - EPS) > 0)
                {
                    // Traverse v.
                    for (int k = start; k < end; k++)
                    {
                        qi[nz] = hi[k];
                        qx[nz] = a * hx[k];

                        nz++;
                    }
                }

                // Store the nz count of current column.
                qp[hi[j] + 1] = nz;
            }

            qp[columns] = nz;

            // There might be gaps in the column indices of Q.
            nz = qp[i];
            for (int j = i + 1; j < columns; j++)
            {
                if (qp[j] > 0)
                {
                    nz = qp[j];
                }
                else
                {
                    qp[j] = nz;
                }
            }

        }

        static bool CheckStorageRequirements(int rows, int cols, int nnz)
        {
            const int NZ_THRESHOLD = 100 * 100; // Matrix dimensions > 100 x 100

            if (nnz > NZ_THRESHOLD && nnz > rows * cols / 2)
            {
                Console.WriteLine("WARNING: matrix Q will be dense.");
            }

            const int MAX_BYTES = 200 * 1024 * 1024; // 200 MB

            return ((cols + 1) * 4 + nnz * 4 + nnz * 8) < MAX_BYTES;

        }

        static void PrintMatrix(string caption, SparseMatrix A)
        {
            Console.WriteLine(caption);

            for (int y = 0; y < A.RowCount; y++)
            {
                for (int x = 0; x < A.ColumnCount; x++)
                {
                    Console.Write("{0,5:0.###}", A.At(y, x));
                }
                Console.WriteLine();
            }

            Console.WriteLine();
        }
    }
    
    static class ExtensionMethods
    {
        public static void CopyTo(this SparseMatrix source, SparseMatrix target)
        {
            int m = source.RowCount;
            int n = source.ColumnCount;
            int nnz = source.NonZerosCount;

            if (target.RowCount != m || target.ColumnCount != n)
            {
                throw new ArgumentException("Invalid dimensions.", nameof(target));
            }

            if (target.NonZerosCount < nnz)
            {
                target.RowIndices = (int[])source.RowIndices.Clone();
                target.Values = (double[])source.Values.Clone();
            }
            else
            {
                Array.Copy(source.RowIndices, 0, target.RowIndices, 0, nnz);
                Array.Copy(source.Values, 0, target.Values, 0, nnz);
            }

            source.ColumnPointers.CopyTo(target.ColumnPointers, 0);
        }
    }
}

Hello, thanks for the code, it's really helpful However, it seems that row or column permutations are not used to construct Q, and I don't understand why? I also tried on a rank deficient matrix where some rows are identical,

1.0  3.0  1.0  0.0
2.0  0.0  2.0  0.0
2.0  0.0  2.0  0.0
2.0  0.0  2.0  0.0

and the CreateQ returns a 5x5 matrix instead of a 4x4 matrix, maybe this has something to do with this?
Thank you again for your work!
Best,
Romain

Actually there's nothing wrong here. CSparse will add fictitious rows if the matrix is structurally rank deficient (you can check this by inspecting the value of m2 in SymbolicFactorization). Permutations are applied after Q is computed.

But keep in mind, that the above code is just for demonstration purpose. I haven't used it in production and wouldn't recommend doing so.