How to make cubic spline interpolation algorithm faster?

Question

I realize that when I use "la_solve" I get the Malloc error due to a limitation in the number of samples "la_solve" can handle.

I noticed that if I use about 900 samples, the program runs well.

Therefore, I modified my code to look at entire data in sequential chunks of 901 samples if samples > 901.

However, I used rather naive techniques and am sure this can be made more efficient and robust: (1) I need a condition that if samples are < 900, do not chunk; (2) take care of repetitions in x and y values as I append chunked data; and (3) make the algorithm more efficient, fast, and robust.

I mostly need help in "// Calling the Re-sampler FUNCTION" section at the end of the code

Thanks !!!

import UIKit
import Accelerate.vecLib.LinearAlgebra

public class Resample {

public class func resample(xi: [Double], xii: [Double], a: [Double]) ->[Double]
{
    // xi - original time stamps (original x)
    // xii - desired time stamps (desired x)
    // a - orignal y values
    // ->[Double] - interpolated y values

    let ni = xii.count // Desired x count

    let n = xi.count // Actual x count

    var h: [Double] = Array(repeating:0, count:n-1)

    for j in 0..<n-1 {
        h[j] = (xi[j+1] - xi[j])
    }

    var A: [Double] = Array(repeating:0, count:(n)*(n))

    A[0] = 1.0

    A[(n*n)-1] = 1.0

    for i in 1..<(n-1) {
        A[(n+1)*i-1] = h[i-1]
        A[(n+1)*i] = 2*(h[i-1] + h[i])
        A[(n+1)*i+1] = h[i]
    }

    var b: [Double] = Array(repeating:0, count:n)

    for i in 1..<n-1 {
        b[i] = (3/h[i])*(a[i+1]-a[i]) - (3/h[i-1])*(a[i]-a[i-1])
    }

    let matA = la_matrix_from_double_buffer(A, la_count_t(n), la_count_t(n), la_count_t(n), la_hint_t(LA_NO_HINT), la_attribute_t(LA_DEFAULT_ATTRIBUTES))

    let vecB = la_matrix_from_double_buffer(b, la_count_t(n), 1, 1, la_hint_t(LA_NO_HINT), la_attribute_t(LA_DEFAULT_ATTRIBUTES))

    let vecCj = la_solve(matA, vecB)

    var cj: [Double] = Array(repeating: 0.0, count: n)

    let status = la_matrix_to_double_buffer(&cj, 1, vecCj)

    if status == la_status_t(LA_SUCCESS) {
        //print(cj.count)
    }
    else {
        print("Failure: \(status)")
    }

    var bj: [Double] = Array(repeating:0, count:n-1)

    for i in 0..<n-1 {
        bj[i] = (1/h[i])*(a[i+1]-a[i]) - (1/3*h[i])*(2*cj[i]+cj[i+1])
    }

    var dj: [Double] = Array(repeating:0, count:n-1)

    for i in 0..<n-1 {
        dj[i] = (1/(3*h[i])) * (cj[i+1]-cj[i])
    }

    var P: [Double] = Array(repeating: 0.0, count: (n-1)*4)

    for i in 0..<n-1 {
        P[(i*4)] = dj[i]
        P[(i*4)+1] = cj[i]
        P[(i*4)+2] = bj[i]
        P[(i*4)+3] = a[i]
    }

    var ai: [Double] = Array(repeating: 0.0, count: ni)


    var jl = Double(1)

    for i in 0..<ni {

        let inter = Double(xii[i])

        while (inter > xi[Int(jl)] && (Int(jl) < n)) {
            jl += 1

        }

        let ind = 4*Int(jl)-4

        var pp = Array(P[ind...ind+3])

        let xx = Double(xii[i]-xi[Int(jl)-1])

        ai[i] = pp[0]*pow(xx, 3) + pp[1]*pow(xx, 2) + pp[2]*xx + pp[3]
    }
    return(ai)
}
}



 // Calling the Re-sampler FUNCTION

let x = Array(stride(from: 0, through: 16649, by: 1.0)) // Orig X (Given, say)

let y = Array(stride(from: 0, through: 16649, by: 1.0)) // Orig Y (Given, say)

let f = 2.0 // Hz (Desired Sampling Freq., say)

let g = Array(stride(from: 0, through: x.count-1, by: 900))
// Helper array for chunking (901 samples each time, say)

var xxx = [Double]() // Results for appended chunks X

var yyy = [Double]() // Results for appended chunks Y

for i in 0..<g.count-1 { // loop through

let xc = Array(x[Int(g[i])...Int(g[i+1])]) // x chunk (count = 901)

let yc = Array(y[Int(g[i])...Int(g[i+1])]) // y chunk (count = 901)

let xci = Array(stride(from: xc[0], through: xc[xc.count-1], by: 1.0/f))
// Desired time stamps for chunk of 901 samples

let yci = Resample.resample(xi: xc, xii: xci, a: yc) // Interpolate chunk

xxx += xci // Append interpolation X

yyy += yci // Append interpolation X

}

if(Int(g[g.count-1])<x.count){
// If helper function last index < original no. of samples

let glb = (Int(g[g.count-1])) // last chunk begin index

let gle = (x.count) // last chunk end index

let xce = Array(x[glb...gle-1]) // last x chunk

let yce = Array(y[glb...gle-1]) // last y chunk

let xcei = Array(stride(from: xce[0], through: xce[xce.count-1], by: 1.0/f))
// Desired time stamps for last chunk of samples

let ycei = Resample.resample(xi: xce, xii: xcei, a: yce) // Interpolate last chunk

xxx += xcei // Append interpolation X for last chunk

yyy += ycei // Append interpolation X for last chunk
}

print(xxx) // would like to identify repeated x values and remove them
print(yyy) // remove y values corresponsding to repated x values (as found above)

// Calling the Re-sampler FUNCTION

Answer 1

Note that you are seeking for every ind using binary search, while it's value is monotonically increasing. So you can just check whether next x-value remains in the same interval or moves to the next one. Pseudo code:

 jl = 1  //before cycle

 //inside cycle
  while (current_x > xi[jl]) and (jl < n)
    jl++
  ind = 4 * jl - 4

Is la_solve general method for linear eq.system solving? If yes, then you'd better use dedicated approach for tridiagonal systems (cubic complexity vs linear complexity). Classic text. Python implementation. Wiki. It is important for long input array.