Transformation matrice vecteur (octave)

Loly

Salut à tous
Dans le cadre d'un tp d'analyse matricielle, je dois créer une fonction sur Octave qui permet, à partir d'une matrice [de] $M_{ m\times p}(\R)$ de transformer celle-ci en un vecteur de taille $m\times p$. Autrement dit, la première ligne de cette matrice deviendra le premier bloc du vecteur et ainsi de suite..
Voici mon code, je ne sais pas quoi mettre dans le V() (ligne 5)
Quelqu'un aurait une idée ?
Merci d'avance.

tatof

Si j'ai bien compris tu pars d'une matrice M et tu la transformes en vecteur colonne. De plus tu veux que la première de la matrice soit le premier bloc du vecteur. Dans ce cas, inutile de faire une fonction avec boucle. Utilise la commande:

M2=M';
V=reshape(M2,size(M,1)*size(M,2),1);

J'utilise une transposée (symbole ') car la commande "reshape" prend les éléments ligne par ligne. Du coup, en transposant, la première ligne devient la première colonne. En appliquant "reshape", cette première colonne devient le premier bloc du vecteur. Pour finir, si tu avais une matrice à élément complexe, il vaut mieux utiliser la commande .' d'Octave car la commande ' donnera la matrice transposée conjuguée.

Loly

Merci pour votre réponse.
Vous avez bien compris ma question mais j'ai manqué de clarté.
En effet ici, la boucle me permet de numéroter ces blocs. Je m'en sers ensuite pour résoudre une équation de type $Ax=b$, par la méthode de relaxation par blocs (en algèbre).
Auriez-vous une idée avec une boucle ?

tatof

Dans mon message précédent, $V$ est un vecteur colonne avec un nombre de ligne égale à $mp$. Tu dis que les lignes de ta matrice $M$ sont les différents blocs du vecteur $V$. C'est-à-dire que $V=[\textbf{v}_1,...,\textbf{v}_m]^T$ avec $\textbf{v}_i$ de longueur $p$.

J'ai fais le code pour trouver les différents blocs $\textbf{v}_i$ sous forme d'une fonction:

function [V,V_bloc]=trans_vect(M,num_bloc)

%% Entrées
% M: matrice de taille mxp
% num_bloc: numéro du bloc à extraire du vecteur V

%% Sorties
% V: vecteur colonne où la première ligne de la matrice M est le premier bloc de V
% V_bloc: bloc extrait du vecteur V, dépend de la variable num_bloc

%% Vérification numéro du bloc
if (num_bloc>size(M,2))
    disp('Erreur dans la sélection, num_bloc<size(M,2)');
end

%% Transformation de M en vecteur colonne
M2=M.';
V=reshape(M2,size(M,1)*size(M,2),1);

%% Extraction d'un bloc
V_bloc=V(1+(num_bloc-1)*size(M,2):num_bloc*size(M,2));

end %fin fonction

Pas besoin de boucle pour la numérotation. Tu utilises juste le nombre de colonne de $M$. Dernière chose, tu parles de relaxation par bloc dans ton dernier message. Es-tu en train de faire la méthode de Gauss-Seidel ou Jacobi? Si oui, j'ai déjà écris un code Octave sur ces 2 méthodes. Tu peux les retrouver dans les archives du forum.

Loly

J'utilise la méthode de Gauss-Seidel, peux-tu m'envoyer le lien direct vers les archives ?

Comment ferais-tu pour une fonction réciproque, c'est dire à partir d'un certain vecteur, le transformer en matrice de taille $m\times p$...? Je sais qu'il est important de vérifier que $m$ est divisible par $p$.

Merci pour ton aide ! ;-)

tatof

Pour commencer voici le lien, c'est le dernier message du sujet. J'ai implémenté la méthode de Gauss-Seidel et la méthode de Jacobi dans une fonction: méthode de Gauss-Seidel

Pour ta question sur la fonction réciproque, tu utilises encore la commande "reshape".

function [M]=trans_mat(V,m,p)

%% Entrées
% V: vecteur colonne
% m,p: nombres entiers

%% Sortie
% M: matrice de taille mxp

%% Vérification dimension
if ((m*p ~=size(V,1)) || (mod(m,p)~=0))
    disp('Erreur dans la sélection de m et p');
end

%% Transformation en matrice
M=reshape(V,m,p).'; %on transpose car le premier bloc du vecteur V est la première ligne de M

end %fin fonction

Loly

Voici la méthode de relaxation "optimisée ici" que je dois programmer, avec l'utilisation du pseudo résidu.
J'ai un peu de mal avec les indices... J'ai commencé par une boucle for, dans laquelle j'ai ajouté des boucles if...
Aurais tu une idée ?
Merci.

tatof

Quelles sont les définitions et les dimensions de tes variables sur l'image?

Loly

G est de taille carrée, ici $m$.
X et B de taille $p$ mais je dois les transformer en matrice $M_{m\times p}(\R)$ par la fonction que tu m'as proposé si dessus.
$w$ et $gamma$ sont en paramètre.

tatof

Pour clarifier: $X=[X_1,...,X_p]^T$ est un vecteur colonne de $ \mathbb{R}^{mp \times 1}$ et chaque $X_i$ est un vecteur colonne de $ \mathbb{R}^{m \times 1}$.

tatof

Le mieux pour commencer est d'écrire quelques termes de ton algorithme pour vérifier si tu le comprends. Cela te permettra de voir également comment mettre les boucles "for" et "if" dans l'algorithme.

Comme hypothèse pour les dimensions de $B$ et $X_{\text{ini}}$, j'ai utilisé mon message précédent en espérant que j'ai bien compris ton problème. De plus, je n'ai pas utilisé les fonctions données dans les précédents messages. Il est inutile d'avoir ces fonctions qui correspondent juste à un "reshape". Voici le code:

function [X_sol]=relax_opti(G,gamma,omega,B,X_ini,p,iter_max)

%% Entrées
% G: matrice réelle de taille mxm vérifiant m divise size(B,1)
% gamma: paramètre réel
% omega: paramètre réel
% B: vecteur colonne de taille mp avec B=[B1,...,Bp]^T
% X_ini: vecteur colonne d'initialisation de taille mp avec X_ini=[X_ini1,...,X_inip]^T
% p: paramètre de découpage vérifiant p divise size(B,1)
% iter_max: nombre d'itération maximale

%% Sortie
% X_sol: vecteur colonne solution de taille mp

%% Vérification taille
if (mod(size(B,1),p)~=0)
    disp('Erreur dans le choix de p');
end

if (mod(size(B,1),size(G,1))~=0)
    disp('Erreur dans les dimensions de G');
end

%% Méthode de relaxation optimisée
X_ini2=reshape(X_ini,size(X_ini,1)/p,p); %mise sous forme de matrice de taille mxp
X_b3=[];
for iter=1:iter_max
    if (iter>2)
        X_b3(1:size(X_ini2,1),:)=[]; %on garde uniquement l'itération précédente
    end
    X_b2=[];
    for i=1:p
        
        %Sélection des blocs
        B_vec=B(1+(i-1)*(size(B,1)/p):i*(size(B,1)/p));
        if (iter==1)
            if (i==p)
                X_b=X_ini2(:,i);
            else
                X_b=X_ini2(:,i:i+1);
            end
        else 
            if (i==p)
                X_b=X_b3(:,i);
            else
                X_b=X_b3(:,i:i+1);
            end
        end
        
        %Pseudo-résidus
        if (i==1)
            R=-G*X_b(:,1)-gamma*X_b(:,2)+B_vec;
        elseif (i==p)
            R=-G*X_b(:,1)-gamma*X_cal+B_vec;
        else
            R=-G*X_b(:,1)-gamma*X_cal-gamma*X_b(:,2)+B_vec;
        end
        
        %% Résolution
        X_cal=X_b(:,1)+omega*inv(G)*R;
        X_b2=[X_b2 X_cal]; %matrice de taille mxp
    end
    X_b3=[X_b3;X_b2];
end
X_sol=reshape(X_b3(end-size(X_ini2,1)+1:end,:),size(B,1),1);
end %fin fonction

Loly

Voici le sujet de mon TP et ce que j'ai fait jusqu'ici. La partie 5 n'étant pas à traiter..
Je ne saisis pas trop ce qu'est ici X_b..
J'ai peur que ton script soit un peu trop complexe, serais tu capable de le simplifier ?

tatof

Pour commencer, je vais parler de ton code. Déjà, il y a quelque chose qui ne va pas dans ta boucle "while". Tu utilises la condition "&&" or cela devrait être la condition ou "||". Sinon tu ne te sers jamais de la condition sur l'erreur. Tu vas toujours aller à ton nombre d'itération maximale. Deuxième chose, personnellement, je n'aime pas trop les boucles "while", je préfère utiliser des boucles "for" car les boucles "while" peuvent tourner indéfiniment si on oublie d'ajouter le compteur. Ici tu as bien mis le compteur sur les itérations. Autre chose, ton calcul d'erreur est sûrement faux puisque tu poses $\text{err}=\text{seuil}+1$. Normalement, cela devrait être la solution trouvée à une itération moins la solution voulue. Tu pourras voir cela dans le code que j'ai joint à ce message. Ton calcul d'erreur doit aussi être dans la boucle sinon il ne sert à rien. Dans mon code j'ai opté pour une boucle "for" sur les itérations et une boucle "break" me permettant de sortir de la boucle des itérations si la condition sur l'erreur est remplie.

A présent par rapport au code que j'ai écris dans mon précédent message. Dans ton code, je vois que tu n'as pas mis les sauvegardes que j'ai fais. Pourtant tu en as besoin car dans le calcul des résidus, tu utilises l'itération précédente pour une colonne fixée. De plus, la variable X_b me permet de choisir entre le vecteur d'initialisation pour la première itération et le vecteur mémoire correspondant à l'itération précédente, ici la variable X_b3.

Pour finir, voici le code complet en utilisant le format donné dans ton tp. J'ai gardé le même nom de fonction que toi "relax_bloc". Je joins également un fichier test correspondant à la partie $4$ de ton tp. Il te permettra d'essayer pour différentes initialisations.

function [X,iter]=relax_bloc(m,p,alpha,beta,gamma,B,omega,X0,seuil,kmax)

%% Entrées
% m: paramètre réel égal à size(B,1)
% p: paramètre réel égal à size(B,2)
% gamma: paramètre réel
% B: matrice bloc de taille mxp avec B=[B1,...,Bp]
% omega: paramètre réel
% X0: matrice bloc d'initialisation de taille mxp avec X0=[X0_1,...,X0_p]
% seuil: seuil fixé pour l'erreur
% kmax: nombre d'itération maximale

%% Sortie
% X: matrice bloc solution de taille mxp
% iter: nombre d'itération réalisée 

%% Matrice G
v1=ones(m,1)*alpha;
v2=ones(m-1,1)*beta;
G=diag(v1)+diag(v2,1)+diag(v2,-1);

%% Méthode de relaxation optimisée
X_b3=[];
for iter=1:kmax
    if (iter>2)
        X_b3(1:size(X0,1),:)=[]; %on garde uniquement l'itération précédente
    end
    X_b2=[];
    for i=1:p
        
        %Sélection des blocs
        if (iter==1)
            if (i==p)
                X_b=X0(:,i);
            else
                X_b=X0(:,i:i+1);
            end
        else 
            if (i==p)
                X_b=X_b3(:,i);
            else
                X_b=X_b3(:,i:i+1);
            end
        end
        
        %Pseudo-résidus
        if (i==1)
            R=-G*X_b(:,1)-gamma*X_b(:,2)+B(:,i);
        elseif (i==p)
            R=-G*X_b(:,1)-gamma*X_cal+B(:,i);
        else
            R=-G*X_b(:,1)-gamma*X_cal-gamma*X_b(:,2)+B(:,i);
        end
        
        %% Résolution
        X_cal=X_b(:,1)+omega*inv(G)*R;
        X_b2=[X_b2 X_cal]; %matrice de taille mxp                
    end
    X_b3=[X_b3;X_b2];
    err=max(abs(reshape(X_b2,m*p,1)-reshape(B,m*p,1))); %formule de l'erreur
    if (err<seuil)
        break; %on sort de la boucle for iter si la condition sur l'erreur est remplie
    end
end
X=X_b3(end-size(X0,1)+1:end,:);
end %fin fonction

Fichier test correspondant à la partie $4$ du tp:

%test 
clear all;close all;clc;

%% Paramètres
a=2; %bord du carré
b=1; %bord du carré
m=3;
p=4;
delta_x=a/(m+1);
delta_y=b/(p+1);
alpha=(2/(delta_x)^2)+(2/(delta_y)^2);
beta=-(1/(delta_x)^2);
gamma=-(1/(delta_y)^2);
B=ones(m,p);
omega=0.1;
X0=zeros(m,p);
seuil=1e-3;
kmax=10;

%%
[U,iter]=relax_bloc(m,p,alpha,beta,gamma,B,omega,X0,seuil,kmax);

%% Tracé
xt=linspace(0,a,m+2);
yt=linspace(0,b,p+2);
[Xt,Yt]=meshgrid(xt,yt);

figure;
Uplot=[zeros(1,m+2);zeros(p,1) U' zeros(p,1);zeros(1,m+2)];
surf(Xt,Yt,Uplot);
shading interp;

Edit: Par rapport au message d'origine: modification de la formule de l'erreur par celle donnée dans la partie validation de ton tp. Simplification du code avec suppression d'une variable inutile.

Rescassol

Bonjour,

Magnifique ! Copié/collé tel quel dans Matlab 2021, ça marche !
Dommage que le TP exclut la validation, ça aurait été intéressant.
Loly, si tu as d'autres sujets de TP dans cette matière, ça m'intéresse.
Qu'est ce qu'il y a d'autre dans ce cours ?

Cordialement,

Rescassol

Loly

Merci à toi tatof, tu m'as été d'une grande aide.

Loly

Salut,
je me permets de te solliciter à nouveau si cela ne te dérange pas...
Je dois implémenter la factorisation QR, pour le coup tout marche, sauf que R et Q sont échangées (i.e Q est triangulaire sup alors que c'est R qui doit l'être....) quand j'effectue un test.
Voila ce que j'ai fait, j'ai déjà optimisé l'algorithme.

tatof

Voici le code que j'ai fait à partir du pseudo-code que tu as donné. Il est identique à celui que tu as fait sauf que chez moi j'ai mis la variable H. J'ai également testé le code sur une matrice carré et non carré et j'obtiens bien R qui est une matrice triangulaire supérieure.
Voici le code de la fonction calculant la décomposition QR :

function [Q,R]=qr_prog(A)

%% Entrée
% A: matrice à décomposer sous forme A=QR

%% Sorties
% Q: matrice orthogonale
% R: matrice triangulaire supérieure

%%
[m,n]=size(A);
Q=eye(m);

for k=1:n
    a=A(k:m,k);
    e1=zeros(m-k+1,1);
    e1(1)=1;
    alpha=norm(a);
    v=a-alpha*e1;
    nu=alpha*(alpha-a(1));
    
    if (nu~=0)
        H=eye(m-k+1)-(v*v'/nu);
        A(k:m,k:n)=H*A(k:m,k:n);
        Q(1:m,k:m)=Q(1:m,k:m)*H;
    end
end
R=A;   
end %fin fonction

Fichier test pour la décomposition QR :

%test prog QR
clear all;close all;clc;
%% Paramètre
A=[12 -51 4;6 167 -68;-4 24 -41]; %matrice carré
% A=[2 3;2 4;1 1];

%% QR
[Q,R]=qr_prog(A);

Loly

Mon erreur se trouvait dans mon test, j'ai effectivement échangé R et Q dans ma fonction^^
Autre question, je souhaite ensuite résoudre le problème des moindres carrés, avec les équations normales j'ai réussi, mais avec la factorisation QR, je dois poser R1, qui prends R et qui lui retire les 0 de la partie inférieure, afin qu'elle soit carrée et que je puisse l'inverser si $\det(A)\neq0$.
Sais tu comment définir cette matrice R1?

tatof

Pour rendre la matrice R carré j'ai écris cette fonction:

function [R1]=mat_carre(R)

%% Entrée
% R: matrice triangulaire supérieure

%% Sortie
% R1: matrice avec les éléments non nuls de R

%%
[m,n]=size(R);
if (m==n)
    R1=R;
else
    diff=abs(m-n);
    R1=R(1:end-diff,:);
end

end %fin fonction

Le fichier test:

%test prog QR
clear all;close all;clc;

%% Paramètre
% A=[12 -51 4;6 167 -68;-4 24 -41]; %matrice carré
A=[2 3;2 4;1 1];

%% QR
[Q,R]=qr_prog(A);

%% Rendre R carré
[R1]=mat_carre(R);

Héhéhé

Tu devrais apprendre à faire une capture écran.