Solving unknown matrics to the power 20

Question

yasmin ismail 2024-7-10，17:19

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2136278-solving-unknown-matrics-to-the-power-20

编辑： Torsten 2024-7-15，20:45

% [P] * [TPM]^n* [R]= CR
disp('Create an array with seven elements in a single row:')
disp('>> P = [1 0 0 0 0 0 0]')
p = [1 0 0 0 0 0 0]
disp('Create an array with seven elements in a single column:')
disp('>> R = [9; 8; 7; 6; 5; 4; 3]')
R = [9; 8; 7; 6; 5; 4; 3]
CR = [0.45]

Solve (TPM)^20=CR/ (p*R)

I would like to find the matrics [TPM] with size 7*7? How to do it ?

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Torsten 2024-7-10，17:32

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2136278-solving-unknown-matrics-to-the-power-20#answer_1483928

编辑：Torsten 2024-7-10，17:44

在 MATLAB Online 中打开

p = [1 0 0 0 0 0 0];
R = [9; 8; 7; 6; 5; 4; 3];
CR = [0.45];
n = 20;
TPM = diag([(CR/(p(1)*R(1)))^(1/n),zeros(1,6)]);
p*TPM^n*R
ans = 0.4500

Hint:

Assume TPM as diagonal: TPM = diag([d(1),...,d(7)]).

Then your relation reads

p(1)*R(1)*d(1)^n + ... + p(7)*R(7)*d(7)^n = CR

30 个评论
显示 28更早的评论隐藏 28更早的评论

Torsten 2024-7-10，19:09

编辑：Torsten 2024-7-10，19:29

在 MATLAB Online 中打开

summation of each row =1 or zero

I only know the definition that all rows must add to 1. That's implemented in the code below.

For your case given MATLAB cannot find a solution. But this doesn't mean that no solution exists.

p = [1 0 0 0 0 0 0];
R = [9; 8; 7; 6; 5; 4; 3];
CR = [0.45];
n = 20;
x0 = eye(7);
obj = @(x)(p*x^n*R-CR)^2;
x = fmincon(obj,x0,[],[],[],[],zeros(7),ones(7),@nonlcon)
Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND =  4.956251e-17.
Local minimum possible. Constraints satisfied.

fmincon stopped because the size of the current step is less than
the value of the step size tolerance and constraints are 
satisfied to within the value of the constraint tolerance.
x = 7x7
    0.4208    0.1777    0.2448    0.1360    0.0014    0.0171    0.0022
    0.0023    0.5024    0.2325    0.1762    0.0555    0.0290    0.0022
    0.0023    0.1365    0.5155    0.2118    0.0858    0.0372    0.0110
    0.0023    0.1443    0.2853    0.3782    0.1210    0.0473    0.0215
    0.0023    0.1498    0.3060    0.2865    0.1715    0.0546    0.0292
    0.0023    0.1522    0.3153    0.3008    0.1579    0.0399    0.0315
    0.0023    0.1512    0.3117    0.2966    0.1568    0.0618    0.0195
<mw-icon class=""></mw-icon>
<mw-icon class=""></mw-icon>
p*x^n*R-CR
ans = 6.1420
sum(x,2)
ans = 7x1
    1.0000
    1.0000
    1.0000
    1.0000
    1.0000
    1.0000
    1.0000
<mw-icon class=""></mw-icon>
<mw-icon class=""></mw-icon>
function [c,ceq] = nonlcon(x);
  ceq(1:7) = sum(x,2)-1;
  c = [];
end

yasmin ismail 2024-7-13，18:15

编辑：Torsten 2024-7-13，18:19

在 MATLAB Online 中打开

@Torsten

yes for above structure but I got different answer, I wrote the code like this:

p = [1 0 0 0 0 0 0];
R = [9; 8; 7; 6; 5; 4; 3];
CR = [0.45];
n = 20;
d = rand(6,1);
du = 1- d;
x0 = diag([d;1]) + diag(du,1)
x0 = 7x7
    0.5368    0.4632         0         0         0         0         0
         0    0.7662    0.2338         0         0         0         0
         0         0    0.2375    0.7625         0         0         0
         0         0         0    0.5983    0.4017         0         0
         0         0         0         0    0.0973    0.9027         0
         0         0         0         0         0    0.0484    0.9516
         0         0         0         0         0         0    1.0000
<mw-icon class=""></mw-icon>
<mw-icon class=""></mw-icon>
obj = @(x)(p*x^n*R-CR)^2;
x = fmincon(obj,x0,[],[],[],[],zeros(7),ones(7),@nonlcon)
Local minimum possible. Constraints satisfied.

fmincon stopped because the size of the current step is less than
the value of the step size tolerance and constraints are 
satisfied to within the value of the constraint tolerance.
x = 7x7
    0.7775    0.1563    0.0179    0.0247    0.0010    0.0123    0.0104
    0.0023    0.8959    0.0457    0.0239    0.0172    0.0137    0.0012
    0.3100    0.4383    0.1072    0.0722    0.0326    0.0221    0.0177
    0.1876    0.6243    0.1296    0.0092    0.0226    0.0145    0.0122
    0.3579    0.3974    0.1403    0.0660    0.0298    0.0007    0.0081
    0.3236    0.2757    0.0803    0.0403    0.0228    0.0002    0.2570
    0.1272    0.6831    0.1459    0.0013    0.0216    0.0093    0.0116
<mw-icon class=""></mw-icon>
<mw-icon class=""></mw-icon>
p*x^n*R-CR
ans = 7.4602
sum(x,2)
ans = 7x1
    1.0000
    1.0000
    1.0000
    1.0000
    1.0000
    1.0000
    1.0000
<mw-icon class=""></mw-icon>
<mw-icon class=""></mw-icon>
function [c,ceq] = nonlcon(x);
ceq(1:7) = sum(x,2)-1;
c = [];
end

Torsten 2024-7-13，18:20

编辑：Torsten 2024-7-14，9:56

在 MATLAB Online 中打开

x0 is only an initial guess. It does not guarantee that the solution from "fmincon" will also have this structure. You will have to rewrite the code solving in only 6 unknowns (the elements on the diagonal except the known 1 at the end) instead of 49 unknowns.

But it seems unlikely that restricting the solution space of possible TPM matrices will yield a solution.

p = [1 0 0 0 0 0 0];
R = [9; 8; 7; 6; 5; 4; 3];
CR = [0.45];
n = 20;
lb = zeros(6,1);
ub = ones(6,1);
obj = @(x)(p*(diag([x;1])+diag(1-x,1))^n*R-CR)^2;
x0 = rand(6,1);
x = fmincon(obj,x0,[],[],[],[],lb,ub)
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.
x = 6x1
    0.1715
    0.1762
    0.1705
    0.1813
    0.1661
    0.1718
<mw-icon class=""></mw-icon>
<mw-icon class=""></mw-icon>
% Reconstructed matrix
X = diag([x;1])+diag(1-x,1)
X = 7x7
    0.1715    0.8285         0         0         0         0         0
         0    0.1762    0.8238         0         0         0         0
         0         0    0.1705    0.8295         0         0         0
         0         0         0    0.1813    0.8187         0         0
         0         0         0         0    0.1661    0.8339         0
         0         0         0         0         0    0.1718    0.8282
         0         0         0         0         0         0    1.0000
<mw-icon class=""></mw-icon>
<mw-icon class=""></mw-icon>
error = obj(x)
error = 6.5025

Torsten 2024-7-14，8:50

编辑：Torsten 2024-7-14，13:30

在 MATLAB Online 中打开

Read again:

You will have to rewrite the code solving in only 6 unknowns (the elements on the diagonal except the known 1 at the end) instead of 49 unknowns.

I included how to reconstruct the matrix from the solution vector in the code above.

Again my advice: Learn MATLAB before you try to solve this problem.

As you can see below, the outcome for p*TPM^20*R is always >= 3, so your value of 0.45 is out of reach. I think the components of your R-vector are too large.

Note that the value 3 as smallest value you can get for p*TPM^20*R corresponds to the value of the objective function from the last code I gave you since (3-0.45)^2 = 6.5025.

p = [1 0 0 0 0 0 0];

R = [9; 8; 7; 6; 5; 4; 3];

n = 20;

pr = linspace(0,1,15);

[pr1,pr2,pr3,pr4,pr5,pr6] = ndgrid(pr);

pr1 = pr1(:);

pr2 = pr2(:);

pr3 = pr3(:);

pr4 = pr4(:);

pr5 = pr5(:);

pr6 = pr6(:);

for i = 1:numel(pr1)

d = [pr1(i);pr2(i);pr3(i);pr4(i);pr5(i);pr6(i)];

outcome(i) = p*(diag([d;1])+diag(1-d,1))^20*R;

end

histogram(outcome,'Normalization','probability')

Torsten 2024-7-14，18:24

编辑：Torsten 2024-7-14，18:27

在 MATLAB Online 中打开

This is the last code version with CR = 3 instead of CR = 0.45. But as you can see from the histogram plot, there seem to be many more matrices with p*TPM^n*R=CR apart from the one that "fmincon" determines.

p = [1 0 0 0 0 0 0];
R = [9; 8; 7; 6; 5; 4; 3];
CR = [3];
n = 20;
lb = zeros(6,1);
ub = ones(6,1);
obj = @(x)(p*(diag([x;1])+diag(1-x,1))^n*R-CR)^2;
x0 = rand(6,1);
x = fmincon(obj,x0,[],[],[],[],lb,ub)
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.
x = 6x1
    0.3503
    0.3396
    0.3170
    0.3228
    0.3792
    0.3400
<mw-icon class=""></mw-icon>
<mw-icon class=""></mw-icon>
% Reconstructed matrix
X = diag([x;1])+diag(1-x,1)
X = 7x7
    0.3503    0.6497         0         0         0         0         0
         0    0.3396    0.6604         0         0         0         0
         0         0    0.3170    0.6830         0         0         0
         0         0         0    0.3228    0.6772         0         0
         0         0         0         0    0.3792    0.6208         0
         0         0         0         0         0    0.3400    0.6600
         0         0         0         0         0         0    1.0000
<mw-icon class=""></mw-icon>
<mw-icon class=""></mw-icon>
%Error
error = obj(x)
error = 8.0323e-08

yasmin ismail 2024-7-15，18:28

@Torsten yes i got it

but can you explain to me what does to the power 2 in this equation related to?

obj = @(x)(p*(diag([x;1])+diag(1-x,1))^n*R-CR)^2

Torsten 2024-7-15，18:42

编辑：Torsten 2024-7-15，20:45

Searching for an x such that f(x) = 0 can be done by minimizing f(x)^2.

Minimizing f(x) doesn't work because "fmincon" would try to make f(x) negative up to -Inf which is not desired.

Instead of f(x)^2, you could take f(x)^n for n being any even integer or any function g(f(x)) where g(y) has its minimum in y=0.

请先登录，再进行评论。

Answer 2

Steven Lord 2024-7-10，17:59

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2136278-solving-unknown-matrics-to-the-power-20#answer_1483948

在 MATLAB Online 中打开

p = [1 0 0 0 0 0 0];
R = [9; 8; 7; 6; 5; 4; 3];
CR = [0.45];
RHS1 = CR./(p*R)
RHS1 = 0.0500

Did you perhaps mean to take the outer product of p and R, not the inner product? That would give you a right-hand side that is 7-by-7, but ...

RHS2 = CR./(R*p)
RHS2 = 7x7
    0.0500       Inf       Inf       Inf       Inf       Inf       Inf
    0.0563       Inf       Inf       Inf       Inf       Inf       Inf
    0.0643       Inf       Inf       Inf       Inf       Inf       Inf
    0.0750       Inf       Inf       Inf       Inf       Inf       Inf
    0.0900       Inf       Inf       Inf       Inf       Inf       Inf
    0.1125       Inf       Inf       Inf       Inf       Inf       Inf
    0.1500       Inf       Inf       Inf       Inf       Inf       Inf
<mw-icon class=""></mw-icon>
<mw-icon class=""></mw-icon>

that's basically garbage. You can raise it to the 1/20th power, but GIGO.

RHS2^(1/20)
ans = 7x7
   NaN   NaN   NaN   NaN   NaN   NaN   NaN
   NaN   NaN   NaN   NaN   NaN   NaN   NaN
   NaN   NaN   NaN   NaN   NaN   NaN   NaN
   NaN   NaN   NaN   NaN   NaN   NaN   NaN
   NaN   NaN   NaN   NaN   NaN   NaN   NaN
   NaN   NaN   NaN   NaN   NaN   NaN   NaN
   NaN   NaN   NaN   NaN   NaN   NaN   NaN
<mw-icon class=""></mw-icon>
<mw-icon class=""></mw-icon>

Raising RHS1 to the (1/20) power or taking the 20th root gives you something that looks less like garbage, but it isn't 7-by-7.

sol1 = RHS1^(1/20)
sol1 = 0.8609
sol2 = nthroot(RHS1, 20)
sol2 = 0.8609

Check by raising the solutions to the 20th power and comparing with RHS1.

format longg
[sol1^20; sol2^20; RHS1]
ans = 3x1
        0.0500000000000001
        0.0500000000000001
                      0.05
<mw-icon class=""></mw-icon>
<mw-icon class=""></mw-icon>

Looks close enough for government work.

What is the actual problem you're trying to solve? I'm guessing this is a homework assignment; if it is show us exactly what you're being instructed to do, not what you think you're being instructed to do.

5 个评论
显示 3更早的评论隐藏 3更早的评论

Steven Lord 2024-7-10，18:29

在 MATLAB Online 中打开

You don't have nearly enough information. Let's take a simpler example: consider a system that has three states and that it starts off in the first state.

initialState = [1; 0; 0];

After two steps, it ends up in the third state always.

finalState = [0; 0; 1];

How did it get from initialState to finalState? There are many different possibilities. One possibility is that your system could have gone from the first state to the second state at the first step then the second state to the third at the second step.

transitionMatrix1 = [0 0 0; 1 0 0; 0 1 1]

transitionMatrix1 = 3x3

0 0 0 1 0 0 0 1 1

transitionMatrix1^2*initialState

ans = 3x1

0 0 1

plot(digraph(transitionMatrix1.'))

Or it could have gone from the first state to the third state at the first step then remained in the third state for the second step.

transitionMatrix2 = [0 0 0; 0 0 0; 1 0 1]

transitionMatrix2 = 3x3

0 0 0 0 0 0 1 0 1

transitionMatrix2^2*initialState

ans = 3x1

0 0 1

plot(digraph(transitionMatrix2.'))

How can you tell whether the transition matrix in this case should be transitionMatrix1 or transitionMatrix2?

Now you want to do this for seven states and twenty steps?

If you had more data the Hidden Markov Model functionality in Statistics and Machine Learning Toolbox might be able to help, but with only initial and final states (which I think are your p and R vectors) I'm not sure hmmestimate will be able to do very much.

Steven Lord 2024-7-10，19:40

Do you assume that TPM is unique? It's not, no more than the solution to this problem posed by Cleve Moler is unique. That means the second of Cleve's Golden Rules of computation applies: "The next hardest things to compute are things that are not unique."

You've already told Torsten one additional requirement, though it doesn't provide enough information to make the solution unique. "If i want to get in the solution matrix summation of each row =0 or 1?" So what other requirements do you have for this problem?

yasmin ismail 2024-7-11，5:58

@Steven Lord I dont what exactly you mean by unique TPM, but what i mean that the condition for my elements will not jump two or three states in one time , they should move one state , if it state 4 shoud move to 3 then 2

and yes the summation for each row in TPM should equal to 1 not zero

请先登录，再进行评论。