使用 Accelerate 从 SparseOpaqueFactorization 获取/提取因式分解

Posted 2023-02-25

【中文标题】使用 Accelerate 从 SparseOpaqueFactorization 获取/提取因式分解

【英文标题】：Get/extract factorisation from SparseOpaqueFactorization using Accelerate

【发布时间】：2021-03-25 08:26:43

【问题描述】：

我正在使用 Apples Swift / Accelerate 框架编写一些线性代数算法。所有的工作和求解的 Ax = b 方程都会产生正确的结果（此代码来自苹果示例）。

我希望能够从中提取 LLT 分解

SparseOpaqueFactorization_Double

对象。但似乎没有任何方法可以提取（打印）分解。有谁知道从 SparseOpaqueFactorization_Double 对象中提取分解矩阵的方法？

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")

【参考方案1】：

好的，经过对苹果 swift 标头的大量调查，我已经解决了这个问题。

有一个名为 Accelerate API 的调用

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

返回这个 SparceOpaqueSubfactor_ 类型。这可用于矩阵乘法以产生“透明”结果（即您可以使用/打印/查看的矩阵）。因此，我将 Cholesky 分解的下三角部分的 SubFactor 乘以 Identity 矩阵以提取因子。工作一种享受！

    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

我编写了一个小函数来生成我在上面的代码中使用的单位矩阵：

    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