Solving unknown matrics to the power 20

% [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 ?

Answer 1

[Torsten]

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

[yasmin ismail]

@Torsten

TPM should be matrix size 7*7 not 0.45

so how to solve this unknown matrix?

[Torsten]

The matrix TPM = diag([(CR/(p(1)*R(1)))^(1/n),zeros(1,6)]) is of size 7x7.

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

TPM = 7x7

0.8609 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0

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

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

Of course there are many more solutions because you have 1 equation with 49 unknowns:

p*TPM^n*R = CR

[yasmin ismail]

@Torsten what does zeros(1,6) mean and how we can get different matrix as solutions?

and If i want to get in the solution matrix summation of each row =0 or 1?

Can i get solution like this based on your answer that i can get different solution?

[Torsten]

what does zeros(1,6) mean

diag([(CR/(p(1)*R(1)))^(1/n),zeros(1,6)]) means that TPM is a diagonal matrix with first diagonal element equal to CR/(p(1)*R(1)))^(1/n) and the six following diagonal elements equal to 0.

how we can get different matrix as solutions

What kind of matrices do you have in mind ? Diagonal matrices are most easily computed because of the relation p(1)*R(1)*d(1)^n + ... + p(7)*R(7)*d(7)^n = CR they have to fulfill. In your case, p(i) = 0 for i>1, thus the matrix I gave you and its negative counterpart are the only diagonal solution matrices for your problem.

[yasmin ismail]

@TorstenI am applying markov chain model and the unkown matrix is the transition propability matrix [TPM] which has in each column probability distributuin which mean that: summation of each row equal 1 or maybe zero without negative value so how to get this answer for an example of answer that i am looking for:

.4 .6 0 0 0 0 0

0 .7 .3 0 0 0 0

0 0 0 0 0 0 1

[Torsten]

Why is this matrix 3x7 instead of 7x7 ? It doesn't fit to the dimensions of p and R.

And if you know the answer, why do you ask ? As I said, there are many matrices that fulfill your requirement. I fear I will again choose one that you don't like.

[yasmin ismail]

@Torsten no this is not answer , i just give an example of how the answer I would like to be

summation of each row =1 or zero

[Torsten]

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]

@Torsten

yes summation for each row in TPM should equal to 1 not zero , its written by mistake, so there is a modification in code you written above ?

[Torsten]

No.

The constraint

ceq(1:7) = sum(x,2)-1;

means that all rows should sum to 1.

The other condition (0 or 1) would have been difficult to implement.

[yasmin ismail]

@Torsten I took the result of X but when i make the multiplication of matrices

P*(x^20)*R= 4.5 or 0.45 whatever the value of CR for example

I didnt get the result of CR , which mean the solution is incorrect

How to fix it?

[Torsten]

I wrote

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

This means that either the initial guess for P as the identity matrix was too bad or that no solution exists.

I don't know of a tailored method that either returns P or "problem unsolvable". So I suggested to use "fmincon" as sledgehammer. But as for all optimization problems, there is no guarantee to converge towards the global minimum.

[yasmin ismail]

@Torsten can you show me how to use "fmincon" as sledgehammer because I dont have experience in using coding I just want to find solution for my equation and I looked for any software to solve but I didnt find , someone suggested that matlab could help

so can you help me to fix it?

[Torsten]

I gave you a code using the general-purpose solver "fmincon" to set up your problem. You didn't get a solution. You can try to vary x0 = eye(7) as initial guess for your matrix. But beyond this, I have no idea how to tackle this specific problem.

Usually in Markov theory, the transition matrices are given and don't need to be determined. Do you try to make a kind of regression of the model ?

[yasmin ismail]

@Torsten what is the meaning of :

x0 = eye(7)

unfortunatly I dont have transition matrics i have only the condition rating within 20 years and initial state,

I didnt make regression analysis how can I make it ?

[Walter Roberson]

x0 = eye(7)

means that the variable x0 is to be initialized to

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

[yasmin ismail]

@Walter Roberson how to let x as following as an example:

.7 .3 0 0

0 .5 .5 0

0 0 .8 .2

0 0 0 1

I want the tranisition matrix to be like this and same time size 7*7

[Torsten]

x0 = [.7    .3     0    0
      0    .5    .5    0
      0    0    .8   .2
      0      0    0    1]

x0 = 4x4

0.7000 0.3000 0 0 0 0.5000 0.5000 0 0 0 0.8000 0.2000 0 0 0 1.0000

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

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

But I suggest you first pass an introductory course to MATLAB free of costs before you continue working on your problem:

https://matlabacademy.mathworks.com/details/matlab-onramp/gettingstarted

[yasmin ismail]

@Torsten no I mean in code you wrote x0 = eye(7)

1 0 0 0 0 0 0

0 1 0 0 0 0 0

0 0 1 0 0 0 0

0 0 0 1 0 0 0

0 0 0 0 1 0 0

0 0 0 0 0 1 0

0 0 0 0 0 0 1

I want x0 to be for example like

0.7000 0.3000 0 0

0 0.5000 0.5000 0

0 0 0.8000 0.2000

0 0 0 1.0000

if p=.7 and q=.3 that mean p+q=1

I want x0 in above code to be p+q=1 in each row this how i want solution of transition matric looks like

[Torsten]

So you want a random matrix with the above structure ?

d = rand(6,1);
du = 1- d;
x0 = diag([d;1]) + diag(du,1)

x0 = 7x7

0.2836 0.7164 0 0 0 0 0 0 0.0907 0.9093 0 0 0 0 0 0 0.4637 0.5363 0 0 0 0 0 0 0.0433 0.9567 0 0 0 0 0 0 0.4659 0.5341 0 0 0 0 0 0 0.0409 0.9591 0 0 0 0 0 0 1.0000

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

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

By the way: This is code. It explicitly sets x0 to the specific matrix.

x0 = [.7    .3     0    0
      0    .5    .5    0
      0    0    .8   .2
      0      0    0    1]

x0 = 4x4

0.7000 0.3000 0 0 0 0.5000 0.5000 0 0 0 0.8000 0.2000 0 0 0 1.0000

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

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

[yasmin ismail]

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

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

[yasmin ismail]

@Torsten in this case you got transition size 6*1 which in correct should be 7*7 because p size1*7 and R size 7*1, so how x = 6x1

0.1904

0.1700

0.1558

0.1730

0.1696

0.1646

[Torsten]

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')

[yasmin ismail]

@Torsten yes I agree with you 3 is the minimum not 0.45, thus CR= [3]

but can you please wrote me the construct code after the last modification because I got confused with the others code written

[Torsten]

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]

@Torsten now it gives the correct answer with error = 8.0323e-08

but the answer should be not less than 3

correct?

[Torsten]

From the histogram plot, the minimum value for p*TPM^20*R is 3 and the maximum value is 9 for the given vectors p and R. So you can only get a solution matrix TPM of this structure if you choose CR >= 3 and CM <= 9 (thus the minimum and maximum components of the vector R).

[yasmin ismail]

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

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]

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.

[yasmin ismail]

@Steven Lord I am applying markov chain model , I dont have the transition probability matrix which is unknown

[Steven Lord]

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

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

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

transitionMatrix1^2*initialState

ans = 3x1

0 0 1

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

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

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

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

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

transitionMatrix2^2*initialState

ans = 3x1

0 0 1

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

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

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.

[yasmin ismail]

@Steven Lord I have the initial state which is

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

and I have a vector of condition rating

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

I have number of years n=20

I have the final condition rating after 20 years = 0.45

I dont have the transition probability matrix, and as you know the equation to predict future condition at n (year) = initial state * (TPM^n) * vactor of condition rating

this is my equation and as you see TPM is unknown and i would like to solve it

and my assumption that the element will move one state from condition to onother what i mean if the element in condition 9 for example it will not jump to 7 directly first move to 8 then 7 this is my assumption and based on it i want to get transition propability matrix

[Steven Lord]

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]

@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