Образец данных

Вот некоторые образцы данных, часто используемые для проверки сплайн аппроксимации с переменными узлами, так называемые титановые тепловые данные. Они регистрируют некоторое свойство титана, измеренное как функция температуры.

[xx,yy] = titanium;
plot(xx,yy,'x');
axis([500 1100 .55 2.25]);
title('The Titanium Heat Data');
hold on

Обратите внимание на довольно резкий пик. Мы используем эти данные для иллюстрации некоторых методов выбора узлов.

Во-первых, мы выбираем несколько точек данных из этих несколько грубых данных. Мы проведем интерполяцию с использованием этого подмножества, а затем сравним результаты с полным набором данных.

pick = [1 5 11 21 27 29 31 33 35 40 45 49];
tau = xx(pick);
y = yy(pick);
plot(tau,y,'ro');
legend({'Full Dataset' 'Subsampled Data'}, 'location','NW');

Общие соображения

Сплайн порядка k с n+k узлы имеет n степени свободы. Поскольку у нас есть 12 сайтов данных, tau(1) < ... < tau(12), для посадки с кубическим сплайном, т.е. сплайном четвертого порядка, требуется последовательность узлов t длины 12 + 4.

Более того, последовательность узлов t должны удовлетворять условиям Шёнберга-Уитни, т.е. должны быть таковы, что i-й узел данных лежит в опоре i-го B-сплайна, т.е.

t(i) < tau(i) < t(i+k) for all i,

при равенстве, допускаемом только в случае узла множественности k.

Одним из способов выбрать последовательность узлов, удовлетворяющую всем этим условиям, является оптимальный узел Gaffney/Powell и Micchelli/Rivlin/Winograd.

Оптимальные узлы

При оптимальной интерполяции сплайна к значениям на площадках

tau(1), ..., tau(n)

скажем, узлы выбираются так, чтобы минимизировать постоянную в стандартной формуле ошибки. В частности, первый и последний узлы данных выбраны в виде k-кратных узлов. Остающееся n-k узлы поставляются optknt.

Вот начало справки от optknt:

ОПТКНТ Оптимальное распределение узлов.

OPTKNT(TAU,K) returns an `optimal' knot sequence for

interpolation at data sites TAU(1), ..., TAU(n) by splines of

order K. TAU must be an increasing sequence, but this is not

checked.

OPTKNT(TAU,K,MAXITER) specifies the number MAXITER of iterations

to be tried, the default being 10.

The interior knots of this knot sequence are the n-K

sign-changes in any absolutely constant function h ~= 0 that

satisfies

integral{ f(x)h(x) : TAU(1) < x < TAU(n) } = 0

for all splines f of order K with knot sequence TAU.

Попытка ОПТКНТ

Мы стараемся использовать optknt для интерполяции на нашем примере, интерполяции кубическими сплайнами к данным

(tau(i), y(i)), for i = 1, ..., n.

k = 4;
osp = spapi( optknt(tau,k), tau,y);

fnplt(osp,'r');
hl = legend({'Full Dataset' 'Subsampled Data' ...
        'Cubic Spline Interpolant Using Optimal knots'}, ...
        'location','NW');
hl.Position = hl.Position-[.14,0,0,0];

Это немного обескураживает!

Здесь, отмеченные звёздами, также являются внутренними оптимальными узлами:

xi = fnbrk(osp,'knots'); 
xi([1:k end+1-(1:k)]) = [];
plot(xi,repmat(1.4, size(xi)),'*');
hl = legend({'Full Dataset' 'Subsampled Data' ...
        'Cubic Spline Interpolant Using Optimal knots' ...
        'Optimal Knots'}, 'location','NW');
hl.Position = hl.Position-[.14,0,0,0];

Что случилось?

Выбор узла для оптимальной интерполяции рассчитан на то, чтобы сделать максимум по всем функциям f соотношения

norm(f - If) / norm(D^k f)

как можно меньше, где числитель является нормой ошибки интерполяции, f - If, и знаменателем является норма k-я производная интерполятора, D^k f. Поскольку наши данные подразумевают, что D^k f является довольно большой, ошибка интерполяции вблизи плоской части данных имеет приемлемый размер для такой «оптимальной» схемы.

Фактически, для этих данных обычный кубический сплайновый интерполятор, предоставленный csapi делает довольно хорошо:

cs = csapi(tau,y);
fnplt(cs,'g',2);
hl = legend({'Full Dataset' 'Subsampled Data' ...
        'Cubic Spline Interpolant Using Optimal knots' ...
        'Optimal Knots' 'Cubic Spline Interpolant Using CSAPI'}, ...
        'location','NW');
hl.Position = hl.Position-[.14,0,0,0];
hold off

Выбор узла для аппроксимации наименьших квадратов

Узлы должны выбираться при аппроксимации наименьших квадратов сплайнами. Один из подходов заключается в использовании узлов с равными интервалами для начала, а затем использования newknt с приближением, полученным для лучшего распределения узлов.

Следующие разделы иллюстрируют эти шаги с полным набором данных по теплу титана.

Аппроксимация наименьших квадратов с помощью однородной узловой последовательности

Мы начинаем с однородной последовательности узлов.

unif = linspace(xx(1), xx(end), 2+fix(length(xx)/4));
sp = spap2(augknt(unif, k), k, xx, yy);
plot(xx,yy,'x');
hold on
fnplt(sp,'r');
axis([500 1100 .55 2.25]);
title('The Titanium Heat Data');
hl = legend({'Full Dataset' ...
        'Least Squares Cubic Spline Using Uniform Knots'}, ...
        'location','NW');
hl.Position = hl.Position-[.14,0,0,0];

Это совсем не удовлетворительно. Поэтому мы используем newknt для сплайновой аппроксимации того же порядка и с тем же количеством отрезков полинома, но разрывы лучше распределены.

Использование NEWKNT для улучшения распределения узлов

spgood = spap2(newknt(sp), k, xx,yy);
fnplt(spgood,'g',1.5);
hl = legend({'Full Dataset' ...
        'Least Squares Cubic Spline Using Uniform Knots' ...
        'Least Squares Cubic Spline Using NEWKNT'}, ...
        'location','NW');
hl.Position = hl.Position-[.14,0,0,0];
hold off

Это довольно хорошо. Между прочим, даже одного внутреннего узла меньше в этом случае не хватило бы.