On 16 Wrz, 03:08, Mariusz Marszałkowski <m...@g...com> wrote:
> On 15 Wrz, 21:57, bartekltg <b...@g...com> wrote:> > Danych jest (było gdy
wysyłałem) 10tys.
>
> > Nie zrozumiales. Macierz jest 10 000 na 16. Ale jej
> > rząd to 13. Maceirz ejst osobliwa, liniowo zalezna.
> > To stad w pierwszej kolejnosci biorą sie zypelnie
> > rozne wektroy dajace ten sam wynkik.
>
> Zrozumiałem ale nie mogłem uwierzyć.
> Dane generowałem funkcją Rand z excela,
> jak to możliwe że dla 10tys losowych wektorów
> układ równań nie jest jednoznaczny?

Ktoś tu niedawno ładnie powiedział, ze pomiedzy osobliwoscia
a nieosobliwoscia jest obszar 'numerycznej osobliwosci':)

Ale tu problem jest bardziej algebraiczny.
Patrzmy na wektorki wierszowe. Są dlugości N i mają
k 'przedziałow' o tej własnosci, ze jest w nim dokladnie
jedna jedynka. Popatrzmy na przestrzen rozpieta przez
wszytkie takie wektory (jest ich mniej niz 10tys, dla N=12,
k=4 81, ogolnie, (N/k)^k ).

Popatrzmy teraz na nastepujecy wektor: w jednym przedziale
(całm) +1, w drugim -1, a w pozostalych 0. Taki wektor nie
lezy w naszej przestrzeni (jest prostopadly do wszytkich innych).

[1;1;1; 0;0;0; -1; -1; -1; 0;0;0;] czy
[ 0;0;0; -1; -1; -1; 0;0;0; 1;1;1] sa prostopadle do kazdego

[ X;Y;V;Z ] gdzie XYVZ za trojwymairowymi wektorami z jedna jedynka.

Wymiar takiej przestrzeni to k-1.

Niezaleznie, ile wektorkow wylosujesz, Twoje pokrycie
ma wymiar co najwyzej N- (k-1).

Skoro taka jest struktura zadania, widac, tak być musi:)
Od biedy mozna sie zastanawiać nad usunieciem
niektorych kolumn (rownowaznie, pewne parametry
stale zero), ale ni ma wielkiej potrzeby jesli godzimy
sie na niejednoznacznosc parametrow.

> > BTW, kopnales sie w wyliczaniach:

> Masz rację, ponieważ nie trzeba mnożyć przez zera.
> Danych jest 10^9 rekordów.
> Rekord ma długość 10^6
> Rekord zawiera np. tylko 4 jedynki i resztę zer.
> Czyli aby ustalić to co w excelu nazwałem
> "poprawką", wystarczy średnio 10^9/10^6*4
> operacji, czyli zaledwie około 4tys.
> "poprawek" jest 10^6, czyli 10^6*4tys...
> Może to się jednak da policzyć dla 4
> niezerowych :)

Oszacowanie na ilosc operacji w celu pojedynczej
poprawki zgadza sie. Ale z tym excelem to dowcip rozumiem;)
Ja wole za podstawową jednostke mieć sprasowany wektor
kolumnowy macierzy A (a nie wierszowy) wiec mam
10^6 rekordow po średnio 10^9/10^6*4 wartosci.
Dzieki temu mam wszytko na widelcu.

Pewnie co jakis czas jednak warto policzyc od nowa Ax-b,
bo poprawki moga nam sie rozjezdzać.

> > Moze poznym wieczorkeim wrzuce kod (od razu ostrzegam,
> > z lenistwa matlabowski;)
>
> Chętnie zobaczę, ale w zupełności wystarczy mi jeśli
> potwierdzisz/zaprzeczysz że mówisz o takiej samej
> metodzie jak w excelu.

Wesja dość robocza, duzo smieci. Wiecej czasu poswieca
na wygenerowanie odpowiedniej postaci danych (w jezykach
blisko bebechow nie bedzie to problemem) niz na same obliczenia.

'Dokłaność' 10^-10, wykonal 25 duzych iteacji
(tu panuje duzy rozrzut) w minute.

Dla N=1000; M=100000; 5 sekund.

Kasujac odpowiedni komentarz dostajesz macierz rzadką
ktora mozna potestowac neizbyt duze zadania.

Pozdrawiam
bartekltg
PS. byly jakies serwisy do umieszczania kodu,
oczywiscie zapomnailem i nie moge znalesc..


function [x,b,y,A ] = main( tol )
%MAIN Summary of this function goes here
% Detailed explanation goes here

N=1000;
M=1000000;
k=4;
N=floor(N/k)*k;

tic
%wygenerowalismy nasza macierz, wiersz to nasz k jedynkowy
%wektor zapisany w postaci sprasowanej: same indeksy,
%gdzie znajduje sie jedynka. Bardziej przyda sie nam inna postać.
%wektor to sprasowany wektor macierza A z zadania.
wA = floor((1+rand(M,k)*(N/k))+ ones(M,k)*diag(0:k-1)*N/k);
toc
tic
wiersz =(1:M)'*ones(1,k);
% kA{j} to sprasowana j-ta _kolumna_ macierzy A
% z zagadnienia min(norm(Ax-b))
%tworzenie kA dalekie od optymalnego, ale dydaktyczne
% a2 to <a*a>, gdzie a to wektor z kA.
A=[];
%A= sparse((1:M)'*ones(1,k),wA,1,M,N); %nie potrzeba, ale porosciej
wygenerowac kA
kA={};
a2 = [];
toc
tic
for j=1:N, % to jest zle, trwa dwa razy dluzej niz wlasciwe
obliczenia.
r=wiersz(wA==j);
kA{j}=r;
a2(j)=length(r);
end;
toc
disp('...for')
tic

clear wA; %szkoda pamieci,

b=rand(M,1); %wektor wynikow.

x=zeros(N,1); %wektorki

y=-b; %w residuum Ax-b
toc

disp('ognia');
tic
popr=1;
cykli=0;
while popr>tol, %powtarzaj, az suma wzglednych poprawek bedzie mala.
%nie do konca dobre kryterium, ale popr >= sum( A'*y ) po wykonaniu
wewn. pętli

perm=randperm(N); % w kazdym duzym obiegu poprawimy kazdy
wspolczynnik.
popr=0;
for j=1:N,
a=perm(j); % a - indeks poprawianego wspolczynnia
v=kA{a}; % sprasowany wektor kolumnowy odpowiadajacy
indeksowi a
delta=sum(y(v))/a2(a); % <y, A(:,a)>/||A(:,a)||^2
x(a)=x(a)-delta;
y(v)=y(v)-delta;
popr=popr+abs(delta/x(a));
end;
popr
cykli=cykli+1;

end;
toc
cykli %a tak dla ciekawosic.