q* = w1 * q1 + w2 * q2 + ... + w2 * qn
normalize q*

start with q* as above
do until convergence
    for each input quaternion i (index)
        diff = q[i] * inverse(q*)
        u[i] = log(diff, base q*)
    //Now perform the linear blend
    adapt := zero quaternion
    weights := 0
    for each input quaternion i
        adapt += weight[i] * u[i]
        weights += weight[i]
    adapt *= 1/weights
    adaptInOriginalSpace = q* ^ adapt    (^ is the power operator)
    q* = adaptInOriginalSpace * q*

count = 1
average_transform = Matrix.Identity(4)
for new_transform in list_of_matrices:
    factor = 1/count
    average_transform = lerp(average_transform, new_transform, factor)
    count += 1