symmetric solutions of linear matrix equations

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.

[Sam Chak]

Oh, you mean something like this?

A   = magic(4);
B   = magic(4);
X   = linsolve(A, B)

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 1.306145e-17.

X = 4x4

1 0 0 0 0 1 0 0 0 0 1 0 0 0 0 1

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

[Dyuman Joshi]

"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]

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

Answer 1

[Bruno Luong]

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

[Paul]

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]

Yes, the problem is

given A, B two (n x n) matrices

minimizing norm(X*A - B, 'fro')

such that X = X.';

Answer 2

[Bruno Luong]

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

% Adjoint of uexpand (E')
function upper = adj_uexpand(tu, X)
X = X + X.';
upper = X(tu);
end

% Model E*A
function Y = symmlprod(upper, tu, A)
Y = uexpand(upper, tu)*A;
end

% Adkoint of symmlprod, A'*E'
function upper = adj_symmlprod(tu, A, Y)
upper = adj_uexpand(tu, Y*A');
end

% Model followed by the adjoint
function res = normalsymmlprod(upper, tu, A)
Y = symmlprod(upper, tu, A);
res = adj_symmlprod(tu, A, Y);
end

Answer 3

[Bruno Luong]

% 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

[Torsten]

X is a symmetric matrix including unkonwn values

Seems to be even harder since X is a mixture of known and unknown values, and I guess the known values should remain unchanged. But it might also be a language problem about what "including" means in this context.

[zeynep ozkayikci]

actually X is similar to X=[X1, X2, X3; vectorize this matrix to solve to XA=b but when convert a matrix in one

X2, X4, X5; column/row vector there are repeating unkown but I do not want them

X3, X5, X6 ]

[Bruno Luong]

Very confusing description

You say earlier X is symmetric and want to solve A*X = B and now X*A = b.

If you want to get effective help you need to clearly describe the problem.

[zeynep ozkayikci]

You are right, sorry for the misleading. I was not aware of the difference. but right one is X*A = b and I want to compute X=b*A^-1

[Sam Chak]

Hi @zeynep ozkayikci,

Could you please update the description in your question above? This will ensure that everyone has a clear understanding of your problem. Alternatively, it would be helpful if you could provide the matrices as a toy problem. This way, we can test if the proposed solution works.

[Bruno Luong]

If you want so "solve" for symmetric X for the matrix equation

X*C = D

in the sense of minimizing the frobenius norm ("least square" norm of vectorized matrix). Transpose the above equation as

C'*X' = D'

meaning

A*Xt = B; with the definitions:

• A := C'
• B := D';
• Xt := X';

The rewritten system has the unknown matrix Xt on the right side of the product. Use my original answer above to solve the rewritten system then transpose the result Xt to get X.

I'm not going to bother to rewrite the algorithm steps specifically for X*C = D, since you keep changing the description of your problem. BTW It's not difficult if you understand the machinery of kronecker, and KKT system used in my answer.

[Bruno Luong]

Code for equation X*A = B

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

2.7194 4.0389 -1.1828 0.1878 5.4956 -1.5598 0.4999 4.1377 -2.4361

<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(A,eye(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 = 3x3

0.4727 1.6572 0.6139 1.6572 0.9317 1.7200 0.6139 1.7200 1.9193

<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.0178

[Torsten]

Maybe reformulating the objective for the 2nd solution can make MATLAB choose "lsqlin" instead of "fmincon" as solver.

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

1.0119 0.6367 -1.0262 1.6773 1.5118 -1.2001 1.5906 -0.0654 -1.2883

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

% Bruno's solution
[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(A,eye(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 = 3x3

0.7902 0.6909 0.8989 0.6909 0.9847 1.2858 0.8989 1.2858 1.0359

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

%2nd solution - less efficient, but easy to set up
prob = optimproblem;
X = optimvar('X',n,n);
prob.Objective = norm(X*A-B,'fro');
prob.Constraints = X == X.';
x0.X = rand(n);
sol = solve(prob,x0)

Solving problem using fmincon.

Local 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.

sol = struct with fields:
    X: [3x3 double]

sol.X

ans = 3x3

0.7902 0.6909 0.8989 0.6909 0.9847 1.2858 0.8989 1.2858 1.0359

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

[Bruno Luong]

@Torsten If we modify the way the objective is written out, it can selects a more appropriate solver

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.1568 -0.2813 0.0320 2.6358 0.3772 -1.1306 0.7032 -0.7556 -0.5532

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

%2nd solution - less efficient, but easy to set up
prob = optimproblem;
X = optimvar('X',n,n);
prob.Objective = sum(reshape(X*A-B,n^2,1).^2,1);
prob.Constraints = X == X.';
x0.X = B/A; %rand(n);
sol = solve(prob,x0)

Solving problem using lsqlin.
The interior-point-convex algorithm does not accept an initial point.
Ignoring X0.

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.

sol = struct with fields:
    X: [3x3 double]

sol.X

ans = 3x3

0.6677 1.1400 1.0304 1.1400 2.3373 0.3333 1.0304 0.3333 0.9890

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

[Torsten]

So easy - I didn't come up with it. Thanks for your help.

Since x0 is ignored, one could remove its use.

[Bruno Luong]

"Since x0 is ignored, one could remove its use."

Oh you are completely right.

Answer 4

[Paul]

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

[Bruno Luong]

". 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.