Suitesparse includes a sparse QR algorithm. The main issue is that (even with
pivoting) the R factor has the same nonzero structure as a Cholesky factor of
A^T A, which is generally much denser than a factor of A, and this degraded
sparsity impacts Q as well.
I wonder if someone would like to contribute a sparse QR to PETSc. It could
have a default implementation via Cholesky QR and the ability to call SPQR from
Suitesparse.
Barry Smith <bsm...@petsc.dev> writes:
> Are the nonzero structures of all the rows related? If they are, one could
> devise a routine to take advantage of this relationship, but if the nonzero
> structures of each row are "randomly" different from all the other rows, then
> it is difficult to see how one can take advantage of the sparsity.
>
>
>
>> On Aug 29, 2023, at 12:50 PM, Thanasis Boutsikakis
>> <thanasis.boutsika...@corintis.com> wrote:
>>
>> Hi all, I have the following code that orthogonalizes a PETSc matrix. The
>> problem is that this implementation requires that the PETSc matrix is dense,
>> otherwise, it fails at bv.SetFromOptions(). Hence the assert in
>> orthogonality().
>>
>> What could I do in order to be able to orthogonalize sparse matrices as
>> well? Could I convert it efficiently? (I tried to no avail)
>>
>> Thanks!
>>
>> """Experimenting with matrix orthogonalization"""
>>
>> import contextlib
>> import sys
>> import time
>> import numpy as np
>> from firedrake import COMM_WORLD
>> from firedrake.petsc import PETSc
>>
>> import slepc4py
>>
>> slepc4py.init(sys.argv)
>> from slepc4py import SLEPc
>>
>> from numpy.testing import assert_array_almost_equal
>>
>> EPSILON_USER = 1e-4
>> EPS = sys.float_info.epsilon
>>
>>
>> def Print(message: str):
>> """Print function that prints only on rank 0 with color
>>
>> Args:
>> message (str): message to be printed
>> """
>> PETSc.Sys.Print(message)
>>
>>
>> def create_petsc_matrix(input_array, sparse=True):
>> """Create a PETSc matrix from an input_array
>>
>> Args:
>> input_array (np array): Input array
>> partition_like (PETSc mat, optional): Petsc matrix. Defaults to None.
>> sparse (bool, optional): Toggle for sparese or dense. Defaults to
>> True.
>>
>> Returns:
>> PETSc mat: PETSc matrix
>> """
>> # Check if input_array is 1D and reshape if necessary
>> assert len(input_array.shape) == 2, "Input array should be 2-dimensional"
>> global_rows, global_cols = input_array.shape
>>
>> size = ((None, global_rows), (global_cols, global_cols))
>>
>> # Create a sparse or dense matrix based on the 'sparse' argument
>> if sparse:
>> matrix = PETSc.Mat().createAIJ(size=size, comm=COMM_WORLD)
>> else:
>> matrix = PETSc.Mat().createDense(size=size, comm=COMM_WORLD)
>> matrix.setUp()
>>
>> local_rows_start, local_rows_end = matrix.getOwnershipRange()
>>
>> for counter, i in enumerate(range(local_rows_start, local_rows_end)):
>> # Calculate the correct row in the array for the current process
>> row_in_array = counter + local_rows_start
>> matrix.setValues(
>> i, range(global_cols), input_array[row_in_array, :], addv=False
>> )
>>
>> # Assembly the matrix to compute the final structure
>> matrix.assemblyBegin()
>> matrix.assemblyEnd()
>>
>> return matrix
>>
>>
>> def orthogonality(A): # sourcery skip: avoid-builtin-shadow
>> """Checking and correcting orthogonality
>>
>> Args:
>> A (PETSc.Mat): Matrix of size [m x k].
>>
>> Returns:
>> PETSc.Mat: Matrix of size [m x k].
>> """
>> # Check if the matrix is dense
>> mat_type = A.getType()
>> assert mat_type in (
>> "seqdense",
>> "mpidense",
>> ), "A must be a dense matrix. SLEPc.BV().createFromMat() requires a
>> dense matrix."
>>
>> m, k = A.getSize()
>>
>> Phi1 = A.getColumnVector(0)
>> Phi2 = A.getColumnVector(k - 1)
>>
>> # Compute dot product using PETSc function
>> dot_product = Phi1.dot(Phi2)
>>
>> if abs(dot_product) > min(EPSILON_USER, EPS * m):
>> Print(" Matrix is not orthogonal")
>>
>> # Type can be CHOL, GS, mro(), SVQB, TSQR, TSQRCHOL
>> _type = SLEPc.BV().OrthogBlockType.GS
>>
>> bv = SLEPc.BV().createFromMat(A)
>> bv.setFromOptions()
>> bv.setOrthogonalization(_type)
>> bv.orthogonalize()
>>
>> A = bv.createMat()
>>
>> Print(" Matrix successfully orthogonalized")
>>
>> # # Assembly the matrix to compute the final structure
>> if not A.assembled:
>> A.assemblyBegin()
>> A.assemblyEnd()
>> else:
>> Print(" Matrix is orthogonal")
>>
>> return A
>>
>>
>> # --------------------------------------------
>> # EXP: Orthogonalization of an mpi PETSc matrix
>> # --------------------------------------------
>>
>> m, k = 11, 7
>> # Generate the random numpy matrices
>> np.random.seed(0) # sets the seed to 0
>> A_np = np.random.randint(low=0, high=6, size=(m, k))
>>
>> A = create_petsc_matrix(A_np, sparse=False)
>>
>> A_orthogonal = orthogonality(A)
>>
>> # --------------------------------------------
>> # TEST: Orthogonalization of a numpy matrix
>> # --------------------------------------------
>> # Generate A_np_orthogonal
>> A_np_orthogonal, _ = np.linalg.qr(A_np)
>>
>> # Get the local values from A_orthogonal
>> local_rows_start, local_rows_end = A_orthogonal.getOwnershipRange()
>> A_orthogonal_local = A_orthogonal.getValues(
>> range(local_rows_start, local_rows_end), range(k)
>> )
>>
>> # Assert the correctness of the multiplication for the local subset
>> assert_array_almost_equal(
>> np.abs(A_orthogonal_local),
>> np.abs(A_np_orthogonal[local_rows_start:local_rows_end, :]),
>> decimal=5,
>> )
