Matlab/Octave i rozklad Cholesky'ego

Witam!
Zwracam się z prośbą o pomoc, muszę napisać program który rozwiązuje regularne liniowe zadanie najmniejszych kwadratow z macierzą A ∈ R [m×n] oraz wektorem prawej strony b ∈ R [m] na dwa sposoby, a więc rozkladem Choleskyego i obrotami Givensa (rzecz jasna bez uzywania funkcji chol i giv). do tej pory probowalem stworzyc kod rozkladu Choleskyego i mam cos takiego:
k=input('Wprowadz parametr k');
A=(2*(rand(2k,k))-ones(2k,k));
B=A'*A;
C=zeros(k,k);
C(1,1)=sqrt(B(1,1));
suma=0;
for i=1:k
for j=1:k
if i==j
suma=C(j,1);
for m=2:(k-1)
C(j,i)=sqrt(B(j,i)-suma^2);
suma=suma+C(j,m);
endfor;
else
suma2=C(j,1)*C(i,1);
for n=2:(k-1)
C(j,i)=(B(j,i)-suma2)/C(i,i);
suma2=suma2+C(j,n)*C(i,n);
endfor;
endif;
endfor;
endfor;
disp(C);
Co tu jest nie tak? Z gory dziekuje!!!
PS wiem ze moje umiejetnosci sa zenujace, jestem poczatkujacy!

Oto błędy, które masz:

1. C(1,1) liczysz dwukrotnie. Raz, przed pętlami, drugi raz wewnątrz pętli dla i=j=1.
2. Pętla wewnętrzna powinna mieć ograniczenie bieżącego indeksu pętli zewnętrznej, a u Ciebie jest to k.
3. Licząc sumę indeksujesz od 2, a nie od 1, oraz do k-1, zamiast do j-1
4. Zastanów się co robisz w pętlach indeksowanych m. Za każdym razem przypisujesz C(j,i) jakąś wartość. I za każdym razem pierwiastkujesz. Dla i==j powinieneś mieć:

      suma=B(j,j);
      for m=1:(j-1)
        suma=suma-C(j,m)^2;
      endfor;
      C(i,j) = sqrt(suma);

Dla i > j wymyśl już sam.

Dzięki za odpowiedź!
Sporo mi się rozjaśniło, poprawiłem kod zgodnie ze wskazówkami ale coś nadal jest nie tak, bo pod diagonalą oblicza same 0, a nad nią nieskonczonosci. Domyślam się, że błąd jest w momencie dzielenia?
Oto co mam teraz, po poprawkach:
k=input('Wprowadz parametr k');
A=(2*(rand(2k,k))-ones(2k,k));
B=A'*A;
C=zeros(k,k);
suma=0;
for i=1:k
for j=1:i
if i==j
suma=B(j,i);
for m=1:(i-1)
suma=suma-C(j,m)^2;
endfor;
C(j,i) = sqrt(suma);
else
suma2=B(j,i);
for n=1:(i-1)
suma2=suma2-C(j,n)*C(i,n);
endfor;
C(j,i)=suma2/C(i,i);
endif;
endfor;
endfor;
disp(C);

el1906 napisał(a):

Dzięki za odpowiedź!
Sporo mi się rozjaśniło, poprawiłem kod zgodnie ze wskazówkami ale coś nadal jest nie tak, bo pod diagonalą oblicza same 0, a nad nią nieskonczonosci. Domyślam się, że błąd jest w momencie dzielenia?
Oto co mam teraz, po poprawkach:
k=input('Wprowadz parametr k');
A=(2*(rand(2k,k))-ones(2k,k));
B=A'*A;
C=zeros(k,k);
suma=0;
for i=1:k
for j=1:i
if i==j
suma=B(j,i);
for m=1:(i-1)
suma=suma-C(j,m)^2;
endfor;
C(j,i) = sqrt(suma);
else
suma2=B(j,i);
for n=1:(i-1)
suma2=suma2-C(j,n)*C(i,n);
endfor;
C(j,i)=suma2/C(i,i);
endif;
endfor;
endfor;
disp(C);

W przypadku i =/=j zmień wszystkie i na j oraz j na i.

Wielkie dzięki, bardzo mi pomogłeś! Program działa!
Teraz próbuję napisać to samo metodą obrotów Givensa. Pierwsza wartość jest liczona dobrze, ale później coś się psuje, więc pewnie mam znowu błąd w pętlach.
G=input('Wprowadz macierz G');
n=rows(G);
T=eye(size(G));
H=G;
for i=2:n
for j=1:i-1
r=sqrt(((H(i,j))^2)*(H(i-1,j)^2));
c=H(i-1,j)/r;
s=((-1)H(i,j))/r;
T(i,j)=s
T(i-1,j)=c
T(i-1,j+1)=-s
T(i,j+1)=c
H=TH
T=eye(size(G))
endfor
endfor
Próbowałem mnożyć macierze też za pętlą ale również nie wychodzi.

Przede wszystkim, nie próbuj uzyskać tego samego co w rozkładzie Choleskiego, bo ta metoda do tego nie służy.

Błędy:

1. W każdej iteracji powinieneś zainicjalizować nowe T, bo w każdej iteracji masz wyzerować tylko 1 element, Ty modyfikujesz T z poprzednich iteracji.
2. Nie rozumiesz, które elementy macierzy H biorą udział w procesie. Jeśli chcesz wyzerować element H(i,j), przy i>j, to do wyliczenia r potrzebujesz H(j,j) i H(i,j)
3. Podobnie, nie rozumiesz, które elementy macierzy T należy wyliczyć. Elementy c idą na przekątną, czyli indeksy powinny być (i,i) i (j,j). Element s idzie w to miejsce, które chcesz wyzerować, a -s w miejsce o zamienionych indeksach. Najlepiej to sobie rozrysuj na kartce.

Nie próbuj zgadywać co trzeba robić. Przeanalizuj najpierw jak wygląda zerowanie jednego elementu, i upewnij się, że to rozumiesz. Potem dopiero wprowadź sobie pętle for, żeby wyzerować wszystkie elementy pod przekątną.

hmm, ale przeciez na koncu pętli wpisałem T=eye, wiec chyba liczy od nowa dla kazdej pętli(?).
z indeksami oczywiście masz rację... poprawiłem ale dalej coś jest nie tak
G=input('Wprowadz macierz G');
b=input('Wprowadz wektor prawej strony b');
n=rows(G);
m=columns(G);
if length(b)==n
H=G;
for i=1:n
for j=i+1:m
T=eye(size(G));
r=sqrt(((H(i,i))^2)+(H(j,j)^2))
c=H(i,i)/r
s=((-1)H(i,j))/r
T(i,j)=s
T(i,i)=c
T(j,i)=-s
T(j,j)=c
H=TH
endfor
endfor
else
disp('Wprowadzono macierz i wektor o nieprawidlowych wymiarach');
endif
patrze na to co krok po kroku robi ten program i wygląda to chyba nieźle ale tylko pierwsze przejście pętli zeruje element a ja nie wiem co jest nie tak, albo program liczy złe wartości c i s albo coś jest nie tak przy mnożeniu macierzy...
eh, jutro deadline a ja jestem w lesie, mam niby jakieś rozwiązania od znajomych i je rozumiem, oni to robią troche inaczej, ale nie chce zrobić tak jak oni bo wszyscy możemy mieć problemy. Tak czy inaczej wielkie dzięki za pomoc!!!

1. Niepotrzebnie zmieniłeś indeksy w pętli, miałeś dobrze, czyli:

for i=2:n
   for j=1:i-1

   endfor
endfor

2. c = H(j,j)/r
3. Czemu stosujesz 2 różne wymiary m i n? To się robi tylko dla macierzy kwadratowych. Najlepiej zrób sobie dodatnio-określoną, tak jak w przykładzie z Choleskim.
4. Z T masz rację, nie zwróciłem uwagi, tyle że lepiej inicjalizować T na początku, bo w ostatniej iteracji inicjalizacja będzie zbędna.

Dziękuję ogromnie! w tym c był pies pogrzebany!! W moim poleceniu jest dla macierz mxn, wydaje mi się że można to wykonać pod warunkiem, że liczba wierszy jest wieksza lub rowna liczbie kolumn, ale wiem że ten program działa tylko dla kwadratowych.

A gdyby zamiast mnożyc macierze po prostu w pętli nadpisywać kolejne wyrazy, to jest to do zrobienia dla macierzy mxn gdzie m>n? Może spróbuję w ten sposób?
Tak czy inaczej jeszcze raz bardzo dziękuję za pomoc, teraz już z górki!