Форум » Программистские штучки » Ку » Ответить
Ку
administrator: В попытке разобраться в причинах появления "неба в клеточку" и колебаний около нуля при стабилизаторе [1 0 0] пересчитал интегралы в матрице q0. Заодно понял откуда они взялись. Забегая вперед - в матрицах q1 и q2 порядок их формирования абсолютно тот же, но вместо сплайновых кусков P1 P2 P3 P4 используются их производные (первые либо вторые). Это оказалось именно то, что Андрей Михалыч рисовал нам квардратиками как в тетрисе, когда мы ходили к нему еще втроем с Пашей. [pre2] function[q]=q0rc(n) %роспись матрицы интеграла произведений от самих сплайнов %нетрудно видеть что начальные и конечные элементы зеркально совпадают %четвертая строка заполняемая в цикле повторяется со сдвигом на позицию по %вертикали %но это еще не все!каждая строка идентична своему столбцу, то есть мтрица %заполняется "уголками" (1,2)=(2,1) (1,3)=(3,1) и т п %всего уникальных интегралов в матрице 10. %миниальный n=6 но при нем цикл не проходит ни разу и матрица без средней %колонки остается. полный состав начинается с n=7. теория в книге %А.М.Волкова "Геолкартирование с ЭВМ" стр 65 q=zeros(n); q(1,1)=1/252; %intgr(p4)^2 dx q(1,2)=43/1680; %p3*p4 dx q(1,3)=1/84; %p2*p4 dx q(1,4)=1/5040; %p1*p4 dx q(2,1)=q(1,2); q(3,1)=q(1,3); %q(4,1)=q(1,4); %это действительно так, но заполнится в цикле q(2,2)=151/630; %p4^2 + p3^2 dx q(2,3)=59/280; %p2*p3 + p3*p4 dx q(2,4)=1/42; %p4*p2 + p3*p1 dx q(2,5)=q(1,4); q(3,2)=q(2,3); %q(4,2)=q(2,4); %это действительно так, но заполнится в цикле %q(5,2)=q(2,5); %это действительно так, но заполнится в цикле q(3,3)=599/1260; %p4^2 + p3^2 +p2^2 dx q(3,4)=397/1680; %p1*p2 + p2*p3 + p3*p4 dx q(3,5)=q(2,4); q(3,6)=q(1,4); %q(4,3)=q(3,4);%это действительно так, но заполнится в цикле %q(5,3)=q(3,5);%это действительно так, но заполнится в цикле %q(6,3)=q(3,6);%это действительно так, но заполнится в цикле for i=4:n-3, q(i,i-3)=1/5040; %аналогичен q(1,4) q(i,i-2)=1/42; %аналогичен q(2,4) q(i,i-1)=397/1680; %аналогичен q(3,4) q(i,i)=151/315; %единственный уникальный интеграл в цикле p4^2 + p3^2 +p2^2 + p1^2 dx q(i,i+1)=397/1680; %аналогичен q(3,4) q(i,i+2)=1/42; %аналогичен q(2,4) q(i,i+3)=1/5040; %аналогичен q(1,4) end, q(n,n)=q(1,1); q(n,n-1)=q(1,2); q(n,n-2)=q(1,3); q(n,n-3)=q(1,4); q(n-1,n)=q(n,n-1); q(n-2,n)=q(n,n-2); q(n-1,n-1)=q(2,2); q(n-1,n-2)=q(2,3); q(n-1,n-3)=q(2,4); q(n-1,n-4)=q(2,5); q(n-2,n-1)=q(n-1,n-2); q(n-2,n-2)=q(3,3); q(n-2,n-3)=q(3,4); q(n-2,n-4)=q(3,5); q(n-2,n-5)=q(3,6); [/pre2]
Ответов - 7
administrator: "продолжаем разговор" (с)Карлсон Неожиданно оказалось, что точно так же, как умножаются друг на друга "одномерные" матрицы Q для получения двухмерных, умножаются друг на друга тензоры. Это действие называется тензорным произведением или произведением Кронекера и в MATLAB для него есть специальная функция kron, точнее даже 2 функции (отдельно для разреженных и неразреженных матриц). С учетом ее существования на земном шаре широко известные в узких кругах функции aq00, aq11, aq22 могут быть переписаны так: [pre2] function[d]=aq00(n,m) %формирование стабилизатора "0" d=kron(q0(n),q0(m)); function[st]=aq11(n,m) %формирование стабилизатора "p" ("p"overhnost) d1=kron(q1(n),q0(m)); d2=kron(q0(n),q1(m)); st=d1+d2; function[st]=aq22(n,m) %формирование стабилизатора "k" ("k"rivizna) d1=kron(q2(n),q0(m)); d2=kron(q0(n),q2(m)); d3=kron(q1(n),q1(m)); st=d1+d2+2*d3; %двойка вытекает из правил дифференцирования функции двух переменных [/pre2] Время счета aq22(68,68) в варианте Андрей Михалыча 682 сек, в варианте Butenop'a 2,47 сек, в этом - 0,4. Большие размерности в варианте А.М. и Butenop'а на рабочем компе не считались (out of memory). Через krom самое большое число, при которм еще не было out of memory, было aq22(435,435). Время счета для этого случая составило 8,9 сек.
Прохожий: Бу-го-га. Вариант Butenop'a для неразреженных матриц по сути точно такой же как в kron. А вот для разреженных используется немного другой алгоритм (который, кстати, надо повнимательнее покурить, пригодится), использующий возможности обращения только к непустым элементам матрицы, что как раз и уменьшает вычислительную нагрузку. В задачах большой размерности проявляется прелесть 64-х битной системы, у мну 1000*1000 считалось ~680 секунд.
administrator: Произведение Кронекера (тензорное произведение) возникает потому, что в любом неодномерном случае наш сплайн тензорный, т.е. сплайн по переменной x не зависит от сплайна по переменной y. В теории результирующий сплайн по плоскости пишется как произведение Кронекера этих двух сплайнов по осям координат. Обратил внимание, что в программах a20, a21, a22 фактически используется тензорное произведение векторов, но оно заменено другой процедурой. Вот фрагмент исходной a20: [pre2] for i=1:r(1) ax=a0(ab,xy(i,1),h(1)); ay=a0(cd,xy(i,2),h(2)); a=ax'*ay; c(i,:)=q00(a,n,m); end [/pre2] Здесь умножение двух векторов приводит к матрице, а потом она превращается в строку. Я же сделал так: [pre2] for i=1:r(1) ax=a0(ab,xy(i,1),h(1)); ay=a0(cd,xy(i,2),h(2)); c(i,:)=kron(ax,ay); end [/pre2] Т.к. в kron реализована специальная техника для разреженных матриц, этот вариант быстрее. Здесь встает еще один вопрос. Почему бы не сделать функцио a0 векторизованной, чтобы не пришлось вызывать ее для каждой новой координаты, а всего лишь один раз для всего массива точек? Сделал так: [pre2] function[a]=a0rc(ab,x,h) %Вычисление строк матрицы А для массива точек наблюдения по одной оси (любой) %может быть задана и одна точка, и любое их количество. RCgoff 27.09.12 %Как и Александром в q00, использовано обращение к элементам матрицы по %порядку. %В данном случае считаются сами сплайн-функции. (не производные) x=x'; %так сделано, т.к. a0rc для нескольких точек вызывается из программ, %где несколько наблюдений указаны в столбик. l=(ab(2)-ab(1))/h; n=floor(l)+3; %число узлов с учетом расширения области %добавляется один слева - Р4 поместится в него при x={ab(1)...h} %и один справа, для P2 при x={ab(2)-h...ab(2)} %плюс три - потому, что разность всегда даст на 1 меньше кол-ва. a=sparse(0); m0=floor((x-ab(1))/h)+2; %m0 - вектор номеров узла, справа от которого расположена точка, % для которой считаем строку. Для этих узлов будет считаться P3. xx=(x-(ab(1)+(m0-2)*h))/h; %xx - вектор приведенных в интервал 0-1 сдвигов от j-го %узла до точки. Александрова формула. siz=length (xx); a(n,siz)=0; %А - транспонированная (Mlab считает элементы по столбцам) m=m0 + [0:n:n*siz-1]; %создали основной индексный вектор ntx=m0<n-1; %логическое условие для 2 индексного вектора... k0=m.*ntx; % 2 индексный вектор с паразитными нулями... fnd=find(k0>0); % получили вектор номеров неотрицательных значений... k=k0(fnd); xxk=xx(fnd); %и создали второй индексный вектор и xx для x<ab(2) - замена if'a! a(m-1)=(1-xx).^3/6; %P4 a(m)=(4-6*xx.^2+3*xx.^3)/6; %P3 a(m+1)=(1+3*xx+3*xx.^2-3*xx.^3)/6; %P2 a(k+2)=xxk.^3/6; %P1. Для точки с координатой x=ab(2) полином P1 %равен 0, и при этом он должен располагаться в узле, который больше %никогда не используется. Так что его можно отбросить. Из-за этого %нюанса P1 считается не всегда - по 2-му индексному вектору k. %У Волкова и Александра вилка была длиннее, но та же по сути. a=a'; %расставлены точки при возведении в степень, т.к. нам %нужны поэлементные вычисления. Иначе ошибка по MPOWER [/pre2] С учетом появившейся векторности a0 ее можно вынести из цикла в a20, тогда вычисление матрицы A запишется так: [pre2] ax=a0rc(ab,xy(:,1),h(1)); ay=a0rc(cd,xy(:,2),h(2)); for i=1:r(1) c(i,:)=kron(ax(i,:),ay(i,:)); end [/pre2]
administrator: Дальнейшие изыскания были направлены в область визуализации. Для начала была упрощена "классическая" программа визуализации, закомментированная в Александровой mnku. Вот исходный текст этой визуализации: [pre2] %Идея визуализации состоит в следующем. Считается матрица А для узлов сетки %и умножается на вектор коэффициентов. Полученный вектор прреобразуется в %матрицу. %задание точек с шагом h/d z1=ab(1):h(1)/d:ab(2); z2=cd(1):h(2)/d:cd(2); nz1=length(z1); nz2=length(z2); pp=zeros(nz1*nz2,2); %pp-массив координат узлов for i=1:nz1 pp(nz2*(i-1)+1:nz2*i,1)=z1(i); pp(nz2*(i-1)+1:nz2*i,2)=z2; end a=a20(ab,cd,pp,h); ff=a*z; %превращение в матрицу fn=zeros(nz2,nz1); tmp=1:nz1*nz2; fn(tmp)=ff(tmp); [/pre2] Вот мой вариант [pre2] %задание точек с шагом h/d - RC edition - без цикла z1=ab(1):h/d:ab(2);z2=cd(1):h/d:cd(2); nz1=length(z1);nz2=length(z2); [X,Y]=meshgrid(z1,z2); tmp=1:nz1*nz2; ppx(tmp)=X(tmp); ppy(tmp)=Y(tmp); pp=[ppx' ppy']; %вычисение вектора значений ff=a20rc(ab,cd,pp,h)*z; %превращение вектора значений в матрицу значений fn=reshape(ff,nz2,nz1); [/pre2]
administrator: Ну а дальше была переписала программа vizualiz_th. Ускорение достигнуто за счет того, что для каждой "подсетки" выбирались 16 соответствующих ей числа из вектора коэффициентов и умножение происходило не на весь большой вектор, а только на 16 чисел. Но и умножались на этот 16-числовой вектор не матрицы, равные по размеру выходной, но заполненные нулями, а то, что у Александра называется "значения сплайнов в 2D без нулей". Вот полный исходный текст. [pre2] function MatVal=vizrc(ab,cd,h,d,z) %Функция счёта матрицы значений для точек сетки карты. %Основная идея состоит в уменьшении количества расчётов. Значения базиса %считаются один раз (за весь расчет). Для каждой подсетки выбираются %отвечающие ей 16 значений из вектора коэффициентов и множатся на базис. %Получившиеся значения в подсетке пишутся на свое место в выходной матрице. %Идея - А.Новых 2010. Доработка - Л.Ядренников 03.10.2012 ab=ab-ab(1);cd=cd-cd(1); %попытка сделать счет начинающийся не с нуля d=ceil(d); %Округление коэфициента сгущения вверх, для корректной работы m=floor((cd(2)-cd(1))/h(2))+3; %количество узлов сетки СПЛАЙНОВ по у ax=a0_grdrc(d); %Значения сплайнов в точках по одной оси. ValSpln=kron(ax,ax); %Значения сплайнов в 2D, без нулей ValSplnLast=kron(ax,ax(1,:)); %Значения сплайнов для последней строки %специфика в том, что для последней строки и столбца надо по одной оси %иметь сплайны, посчитанные по коэф. сгущения, а по второй оси - только %первый (в точке, совпадающей с узлом). Для столбца можно взять первые d %строк из основной матрицы, а для строки проще пересчитать %Считаем поклеточно выходную матрицу. %1.задание точек с шагом h z1=(ab(1):h(1):ab(2))/h(1);z2=(cd(1):h(2):cd(2))/h(2); nz1=length(z1);nz2=length(z2); [X,Y]=meshgrid(z1,z2); %2.Цикл счета значений в клетках (подсетках по Александру) MatVal=zeros(nz2*d-(d-1),nz1*d-(d-1)); for i=1:nz1*nz2 %столько итераций, сколько клеток. if and(X(i)<z1(end),Y(i)<z2(end)) %для полных клеток indx=1+kron(ones(1,4),(m*X(i)+Y(i):m*X(i)+Y(i)+3))+kron([0:m:4*m-1],ones(1,4)); %индексный вектор vkpart=z(indx); %часть вектора коэффициентов для текущей клетки - 16 цифр VectVal=ValSpln*vkpart; %рассчитываем значения по всей клетке - вектор MatVal((Y(i)*d)+1:(Y(i)*d+1)+d-1,(X(i)*d)+1:(X(i)*d+1)+d-1)=reshape(VectVal,d,d); %записали его в матрицу клеткой else if Y(i)<z2(end)%для последней колонки indx=1+kron(ones(1,3),(m*X(i)+Y(i):m*X(i)+Y(i)+3))+kron([0:m:3*m-1],ones(1,4));%для последней колонки вектор на 4 короче... vkpart=[z(indx);zeros(4,1)];%...потому что последние 4 коэффициента нули MatVal(Y(i)*d+1:Y(i)*d+d,end)=ValSpln(1:d,:)*vkpart; %пишем в выходную матрицу else if X(i)<z1(end)%для последней строки indx=1+kron(ones(1,4),(m*X(i)+Y(i):m*X(i)+Y(i)+3))+kron([0:m:4*m-1],ones(1,4)); %для последней строки нулями должны быть каждые 4-е элементы в векторе коэфф. indx2=(kron(ones(1,4),[ones(1,3) 0]))'; vkpart1=z(indx(1:end-1)); vkpart=[vkpart1.*indx2(1:end-1);0];%...сделали это MatVal(end,X(i)*d+1:X(i)*d+d)=(ValSplnLast*vkpart)'; %пишем в выходную матрицу. else %для последней точки indx=1+kron(ones(1,3),(m*X(i)+Y(i):m*X(i)+Y(i)+3))+kron([0:m:3*m-1],ones(1,4)); %вектор на 5 короче indx=indx(1:end-1); vkpart=[z(indx);zeros(5,1)]; MatVal(end)=ValSpln(1,:)*vkpart; end %конец ветки последней строки-точки end %конец веток последних строк-столбцов-точек end %конец веток по if вообще end %конец цикла end [/pre2] Ускорение существенно. Например, mnku([0 190],[0 190],[5 5 1;5 6 -6; 7 8 10),0,0,[0 1 0],0,1,1,3 с визуализацией Алескандра считает 281с, из них 206с на визуализацию, а с этой 85с и 5с на визуализацию. С ускоренным визуализатором, с ускоренным вычислителем двумерных матриц Q самый медленный элемент в решении задач большого размера - это решение системы и отрисовка контуров (при большом сгущении). При большом количестве наблюдений время тратится и на обработку наблюдений, составление матрицы А. Все-таки там цикл пока что.
administrator: Даже в одномерном случае от тензоров никуда не уйти. Сплайновый базис - не является ортогональным. Для не ортогональных базисов скалярное произведение задается метрическим тензором. Матрица Q суть метрический тензор пространства сплайнов. Не случайно Андрей Михалыч считал нормы как fQf', а скалярные произведения как fQg' Одновременно это билинейная форма и одновременно это матрица Грама. Собственно, A'A это тоже матрица Грама для матрицы А, про это я понял еще в 10 году. А вот что метричский тензор только сейчас. В этой связи мне еще интереснее связь всего этого дела и регуляризации Тихонова. Ведь когда (для обычных задач) матрица Тихонова - это единичная матрица, одновременно это метрический тензор Евклидова пространства. P.S.че-то не Сидоров, никто на это не обращали внимания P.P.S.метрические тензоры это уже математический аппарат общей теории относительности, забавно
administrator: С позиций дифференциальной геометрии поверхностей роль ку-матрицы и сымсл регуляризации, стабилизаторов получает простое истолкование. Строя карту, мы строим поверхность. Геометрия поверхностей - хорошо разработанный раздел математики. Так вот, выше было написано, что матрица Q суть метрический тензор пространства сплайнов. В теории поверхностей метрический тензор называется "первой квадратичной формой" и далее процитирую Википедию: Первая квадратичная форма или метрический тензор поверхности ― квадратичная форма от дифференциалов координат на поверхности, которая определяет внутреннюю геометрию поверхности в окрестности данной точки. Знание первой квадратичной формы достаточно для вычисления длин дуг, углов между кривыми, площади областей на поверхности. Замечу, что Q - именно квадратичная форма от дифференциалов координат. Таким образом, вид метрического тензора задает все характеристики поверхности. В реальной задаче тип стабилизатора определяет тип поверхности, но мы к нему прибавляем матрицу А, и итог - тоже метрический тензор своего рода, только компромисс между идеалом и реальностью. И теперь я понимаю смысл регуляризации. В предыдущем посте написано: "Ведь когда (для обычных задач) матрица Тихонова - это единичная матрица, одновременно это метрический тензор Евклидова пространства". В общем, регуляризация - это метризация решения. Закономерность, которая вводится в решение через стабилизатор, наделяет решение (многообразие) определенными геометрическими (внутренними геометрическими)свойствами. О как! P.S. на это меня натолкнула кандидатская Плавника, где он и Сидоров показывают, что "небо в клеточку" - решение уравнения Клейна-Гордона. Дале по ассоциативному ряду. P.P.S.че-то не Сидоров, никто на это не обращали внимания (на приложение дифференциальной геометрии к этому делу). Даже, кажется, сплайновики из Новосибирска.
полная версия страницы