Tomm P
Tomm P

Reputation: 825

Get/extract factorisation from SparseOpaqueFactorization using Accelerate

I am writing some Linear Algebra algorithms using Apples Swift / Accelerate framework. All works and the solved Ax = b equations produce the right results (this code is from the apple examples).

I would like to be able to extract the LLT factorisation from the

SparseOpaqueFactorization_Double

object. But there doesn't seem to be any way to extract (to print) the factorisation. Does anyone know of a way of extracting the factorised matrix from the SparseOpaqueFactorization_Double object?

import Foundation
import Accelerate

print("Hello, World!")

// Example of a symmetric sparse matrix, empty cells represent zeros.

var rowIndices: [Int32] = [0, 1, 3,   // Column 0
                           1, 2, 3,   // Column 1
                            2,        // col 2
                            3]        // Col 3
 

// note that the Matrix representation is the upper triangular
// here. Since the matrix is symmetric, no need to store the lower
// triangular.

 var values: [Double]  = [10.0, 1.0 ,      2.5,          // Column 0
                            12.0, -0.3, 1.1,        // Column 1
                                  9.5,              // Col 2
                                       6.0 ]        // Column 3

var columnStarts = [0,      // Column 0
                    3,      // Column 1
                    6, 7,   // Column 2
                    8]      // col 3

var attributes = SparseAttributes_t()
attributes.triangle = SparseLowerTriangle
attributes.kind = SparseSymmetric

let structure = SparseMatrixStructure(rowCount: 4,
                                      columnCount: 4,
                                      columnStarts: &columnStarts,
                                      rowIndices: &rowIndices,
                                      attributes: attributes,
                                      blockSize: 1)


let llt: SparseOpaqueFactorization_Double = values.withUnsafeMutableBufferPointer { valuesPtr in
    let a = SparseMatrix_Double(
        structure: structure,
        data: valuesPtr.baseAddress!
    )

    return SparseFactor(SparseFactorizationCholesky, a)
}


var bValues = [ 2.20, 2.85, 2.79, 2.87 ]
var xValues = [ 0.00, 0.00, 0.00, 0.00 ]

bValues.withUnsafeMutableBufferPointer { bPtr in
    xValues.withUnsafeMutableBufferPointer { xPtr in
        
        let b = DenseVector_Double(
            count: 4,
            data: bPtr.baseAddress!
        )
        let x = DenseVector_Double(
            count: 4,
            data: xPtr.baseAddress!
        )

        SparseSolve(llt, b, x)
    }
}
for val in xValues {
    print("x = " +  String(format: "%.2f", val), terminator: " ")
    
}
print("")
print("Success")

Upvotes: 2

Views: 97

Answers (1)

Tomm P
Tomm P

Reputation: 825

OK so after much sleuthing around the apple swift headers, I have solved this problem.

There is an Accelerate API call called

public func SparseCreateSubfactor(_ subfactor: SparseSubfactor_t, _ Factor: SparseOpaqueFactorization_Double) -> SparseOpaqueSubfactor_Double

which returns this SparceOpaqueSubfactor_ type. This can be used in a matrix multiplication to produce a "transparent" result (i.e. a matrix you can use/print/see). So I multiplied the SubFactor for the Lower triangular part of the Cholesky factorisation by the Identity matrix to extract the factors. Works a treat!

    let subfactors = SparseCreateSubfactor(SparseSubfactorL, llt)
    
    var identValues = generateIdentity(n)
    ppm(identValues)
    
    let sparseAs = SparseAttributes_t(transpose: false,
                                      triangle: SparseUpperTriangle,
                                      kind: SparseOrdinary,
                                      _reserved: 0,
                                      _allocatedBySparse: false)
    
    let identity_m = DenseMatrix_Double(rowCount: Int32(n),
                                        columnCount: Int32(n),
                                        columnStride: Int32(n),
                                        attributes: sparseAs,
                                        data: &identValues)
    
    
    SparseMultiply(subfactors, identity_m) // Output is in identity_m after the call

I wrote a small function to generate an identity matrix which I've used in the code above:

    func generateIdentity(_ dimension: Int) -> [Double] {
    
    var iden = Array<Double>()
    
    for i in 0...dimension - 1 {
        for j in 0...dimension - 1 {
            if i == j {
                iden.append(1.0)
                
            } else {
                iden.append(0.0)
            }
            
        }
        
        
    }

    return iden
}

Upvotes: 3

Related Questions