# symmetric solutions of linear matrix equations

5 visualizaciones (últimos 30 días)
zeynep ozkayikci el 6 de Abr. de 2024
Editada: Bruno Luong el 8 de Abr. de 2024
I have X symmetric matrix , X is similar to X=[X1, X2, X3; X2, X4, X5; X3, X5, X6 ] size is changeable. and want to compute the unkown values so I vectorize the X.
I vectorize X matrix to solve XA=b but when converting a matrix in one column/row vector there are repeating unkowns but I do not want them in my solution matrix. How can I solve it?
Vec(X)*A=b
Vec(X)=b*A^-1
I want to write a matlab function to solve the problem.
##### 3 comentariosMostrar 1 comentario más antiguoOcultar 1 comentario más antiguo
Dyuman Joshi el 6 de Abr. de 2024
"X is a symmetric matrix including unkonwn values"
How many values are unknown?
What are the values in A and B?
"I have A and B matrices. I need to solve X
(AX = B )"
(Assuming compatible dimensions) If X is a vector, A*X will be a vector, which will not be equal to be B, a matrix.
Please share the code you have written yet and the values for A and B.
zeynep ozkayikci el 6 de Abr. de 2024
I want to write a matlab function to solve the problem. It will take 2 row vectors x and y for example and it will solve for this equation: vec(theta)*(kronecker product of x minus kronecker product of y) where theta is a symmetric square matrix. But when I turn theta matrix to the vector form some of the unkown values of theta are repeating since its symmetric. How can I solve this problem

Iniciar sesión para comentar.

Bruno Luong el 7 de Abr. de 2024
Editada: Bruno Luong el 7 de Abr. de 2024
If you have tthe optimization toolbox
% Generate test matrices
n = 3;
A = randn(n);
X=rand(n); X = X+X.';
B = X*A;
% Add small noise to B
B = B + 1e-1*randn(size(B))
B = 3x3
-0.3105 -0.5486 -0.2153 -1.4848 -1.0955 -0.9003 -0.3543 0.1766 0.8062
<mw-icon class=""></mw-icon>
<mw-icon class=""></mw-icon>
ij = nchoosek(1:n,2)-1;
m = size(ij,1);
r = (1:m)';
C = zeros([m n^2]);
C(r + ij * [n; 1] * m) = 1;
C(r + ij * [1; n] * m) = -1;
M = kron(A',eye(n));
X = lsqlin(M,B(:),[],[],C,zeros(m,1));
Minimum found that satisfies the constraints. Optimization completed because the objective function is non-decreasing in feasible directions, to within the value of the optimality tolerance, and constraints are satisfied to within the value of the constraint tolerance.
X = reshape(X,n,n)
X = 3x3
0.1068 0.9432 1.0007 0.9432 1.6786 0.5287 1.0007 0.5287 1.4214
<mw-icon class=""></mw-icon>
<mw-icon class=""></mw-icon>
XA = X*A
XA = 3x3
-0.3045 -0.6087 -0.1258 -1.5224 -1.0886 -0.8758 -0.2886 0.1625 0.7665
<mw-icon class=""></mw-icon>
<mw-icon class=""></mw-icon>
norm(XA-B,'fro')/norm(B,'fro') % error
ans = 0.0602
##### 2 comentariosMostrar NingunoOcultar Ninguno
Paul el 7 de Abr. de 2024
Editada: Paul el 7 de Abr. de 2024
Hi Bruno,
Can you explain the mathematical problem that this code is solving?
It looks like the problem is:
Given A, B each n x n, solve for symmetric X s.t that
B - X*A
is minimized in some sense.
Is that correct?
Bruno Luong el 7 de Abr. de 2024
Editada: Bruno Luong el 7 de Abr. de 2024
Yes, the problem is
given A, B two (n x n) matrices
minimizing norm(X*A - B, 'fro')
such that X = X.';

Iniciar sesión para comentar.

### Más respuestas (3)

Bruno Luong el 8 de Abr. de 2024
Editada: Bruno Luong el 8 de Abr. de 2024
This is a method that use the small "original" linear system. I use pcg since it a linear solver that can accept function handle instead of matrix coefficients.
No toolbox is required. Note the convergence of pcg seems to be fragile without preconditioning.
% Generate test matrices
n = 5;
A = randn(n);
X=rand(n); X = X+X.';
B = X*A;
% Add small noise to B
B = B + 1e-1*randn(size(B));
% Result from expanded kron system
ij = nchoosek(0:n-1,2);
m = size(ij,1);
r = (1:m)';
C = sparse(r, ij * [n; 1] +1, 1, m, n^2) - ...
sparse(r, ij * [1; n] +1, 1, m, n^2);
M = kron(A, speye(n));
Y=[M*B(:); zeros(m,1)]';
M=[M*M', C';
C, zeros(m)];
X=Y/M;
X = reshape(X(1:n^2),n,n) % symmetric
X = 5x5
15.7701 -10.4191 8.7905 9.7548 -3.7987 -10.4191 9.3762 -5.2915 -4.8205 4.6536 8.7905 -5.2915 5.2898 5.1365 -2.2491 9.7548 -4.8205 5.1365 4.4069 -1.7434 -3.7987 4.6536 -2.2491 -1.7434 3.1308
<mw-icon class=""></mw-icon>
<mw-icon class=""></mw-icon>
norm(X - X', 'fro');
XA = X*A; % should be close to B
norm(XA-B,'fro')/norm(B,'fro') % error
ans = 0.0294
% New method starts here tor
% Solve
% X = E*Upper, E is a linear operator that expands upper to symmetric matrix
% X*A = B, minimize "fro" norm sense
tu=triu(true(n)); % only UPPER coefficients of that positions is considered
AtB = adj_symmlprod(tu, A, B); % multiply RHS by (E*A)'
upper = pcg(@(upper) normalsymmlprod(upper, tu, A), AtB);
pcg converged at iteration 15 to a solution with relative residual 4.2e-08.
XX = uexpand(upper, tu)
XX = 5x5
15.7701 -10.4191 8.7905 9.7548 -3.7987 -10.4191 9.3762 -5.2915 -4.8205 4.6536 8.7905 -5.2915 5.2898 5.1365 -2.2491 9.7548 -4.8205 5.1365 4.4069 -1.7434 -3.7987 4.6536 -2.2491 -1.7434 3.1308
<mw-icon class=""></mw-icon>
<mw-icon class=""></mw-icon>
norm(XX*A-B,'fro')/norm(B,'fro') % error
ans = 0.0294
% E operator: Expand triangular part to symmetric matrix
% Note the diagonal is double
function X = uexpand(upper, tu)
X = zeros(size(tu));
X(tu) = upper;
X = X + X.';
end
X = X + X.';
upper = X(tu);
end
% Model E*A
function Y = symmlprod(upper, tu, A)
Y = uexpand(upper, tu)*A;
end
function upper = adj_symmlprod(tu, A, Y)
end
% Model followed by the adjoint
function res = normalsymmlprod(upper, tu, A)
Y = symmlprod(upper, tu, A);
end
##### 0 comentariosMostrar -2 comentarios más antiguosOcultar -2 comentarios más antiguos

Iniciar sesión para comentar.

Bruno Luong el 6 de Abr. de 2024
Editada: Bruno Luong el 6 de Abr. de 2024
% Generate test matrices
n = 5;
A = randn(n);
X=rand(n); X = X+X.';
B = A*X;
% Add small noise to B
B = B + 1e-1*randn(size(B))
B = 5x5
-1.0496 -1.2357 -2.2978 -0.0452 -3.8092 -1.8387 -2.2418 -2.9166 -2.9147 -1.3836 1.2053 1.3404 1.6730 0.3741 2.5606 0.2147 -0.4751 0.2241 0.1728 -0.7981 0.1005 0.4824 0.7771 0.2816 1.1183
<mw-icon class=""></mw-icon>
<mw-icon class=""></mw-icon>
[i,j] = find(triu(ones(n),1));
k = sub2ind([n,n],i,j);
l = sub2ind([n,n],j,i);
m = numel(i);
r = (1:m)';
s = [m n^2];
C=accumarray([r k(:)],1,s)-accumarray([r l(:)],1,s);
M = kron(eye(n),A);
Y=[M'*B(:); zeros(m,1)];
M=[M'*M, C';
C zeros(m)];
X=M\Y;
X = reshape(X(1:n^2),n,n) % symmetric
X = 5x5
1.1044 0.8890 0.9223 0.6105 0.4487 0.8890 1.8126 0.8434 1.2570 1.0306 0.9223 0.8434 1.5951 0.7616 1.1951 0.6105 1.2570 0.7616 0.5013 1.6519 0.4487 1.0306 1.1951 1.6519 1.0718
<mw-icon class=""></mw-icon>
<mw-icon class=""></mw-icon>
norm(X - X', 'fro')
ans = 2.3915e-15
AX = A*X % should be close to B
AX = 5x5
-1.0948 -1.2704 -2.2409 -0.0097 -3.8151 -1.8318 -2.3029 -2.8219 -2.8069 -1.4021 1.2111 1.3671 1.8268 0.2867 2.6107 0.1661 -0.4451 0.2532 0.1078 -0.7411 0.0627 0.4175 0.8071 0.4258 1.1224
<mw-icon class=""></mw-icon>
<mw-icon class=""></mw-icon>
norm(AX-B,'fro')/norm(B,'fro') % error
ans = 0.0402
% It should be better than solving original matrix then symmetrizeing it
XX = A\B;
XX = 1/2*(XX+XX');
norm(A*XX-B,'fro')/norm(B,'fro')
ans = 0.1193
##### 11 comentariosMostrar 9 comentarios más antiguosOcultar 9 comentarios más antiguos
Torsten el 7 de Abr. de 2024
Editada: Torsten el 7 de Abr. de 2024
So easy - I didn't come up with it. Thanks for your help.
Since x0 is ignored, one could remove its use.
Bruno Luong el 7 de Abr. de 2024
"Since x0 is ignored, one could remove its use."
Oh you are completely right.

Iniciar sesión para comentar.

Paul el 7 de Abr. de 2024
Movida: Bruno Luong el 7 de Abr. de 2024
Here's an alternative using mldivide, \ that arrives at the same result for this particular problem. Is it effectively the same solution as yours using lsqlin?
rng(100)
% Generate test matrices
n = 3;
A = randn(n);
X=rand(n); X = X+X.';
B = X*A;
% Add small noise to B
B = B + 1e-1*randn(size(B));
ij = nchoosek(1:n,2)-1;
m = size(ij,1);
r = (1:m)';
C = zeros([m n^2]);
C(r + ij * [n; 1] * m) = 1;
C(r + ij * [1; n] * m) = -1;
M = kron(A',eye(n));
X = lsqlin(M,B(:),[],[],C,zeros(m,1));
Minimum found that satisfies the constraints. Optimization completed because the objective function is non-decreasing in feasible directions, to within the value of the optimality tolerance, and constraints are satisfied to within the value of the constraint tolerance.
X = reshape(X,n,n)
X = 3x3
1.3098 0.0872 1.5716 0.0872 0.4135 1.4363 1.5716 1.4363 0.8750
<mw-icon class=""></mw-icon>
<mw-icon class=""></mw-icon>
XA = X*A;
norm(XA-B,'fro')/norm(B,'fro') % error
ans = 0.0604
XX = (kron(eye(3),A')*insmat(3)) \ reshape(B',[],1)
XX = 6x1
1.3098 0.0872 1.5716 0.4135 1.4363 0.8750
<mw-icon class=""></mw-icon>
<mw-icon class=""></mw-icon>
XX = reshape(insmat(3)*XX,3,3)
XX = 3x3
1.3098 0.0872 1.5716 0.0872 0.4135 1.4363 1.5716 1.4363 0.8750
<mw-icon class=""></mw-icon>
<mw-icon class=""></mw-icon>
function Is = insmat(s)
% reference: Vetter, W.J., "Vector Structures and Solutions of Linear
% Matrix Equations," Linear Algebra and Its Applications, 10, 181-188,
% 1975
% first form the index matrix m
m = tril(ones(s,s));
m(m>0) = 1:(s*(s+1)/2);
m = m + tril(m,-1)';
I = eye(s*(s+1)/2,s*(s+1)/2);
Is = I(m(:),:);
end
##### 1 comentarioMostrar -1 comentarios más antiguosOcultar -1 comentarios más antiguos
Bruno Luong el 7 de Abr. de 2024
Movida: Bruno Luong el 7 de Abr. de 2024
". Is it effectively the same solution as yours using lsqlin?"
Not necessary your method parametrizes the subspace of matrices that meet the constraints X = X.', then solve the leastsqure in term of parametrization.
This can also be done by find the (orthogonal) basis null C, which can be carried out by QR decomosition on C'.
Other alternative is forming a KKT system and solve it.
Not sure lsqlin uses which method.
The KKT method is used by my my second answer.

Iniciar sesión para comentar.

### Categorías

Más información sobre Linear Least Squares en Help Center y File Exchange.

R2019b

### Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by