Solving unknown matrics to the power 20
Show older comments
% [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 ?
Answers (2)
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
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 Comments
yasmin ismail
on 10 Jul 2024
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)])
Of course there are many more solutions because you have 1 equation with 49 unknowns:
p*TPM^n*R = CR
yasmin ismail
on 10 Jul 2024
Edited: yasmin ismail
on 10 Jul 2024
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
on 10 Jul 2024
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
on 10 Jul 2024
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)
p*x^n*R-CR
sum(x,2)
function [c,ceq] = nonlcon(x);
ceq(1:7) = sum(x,2)-1;
c = [];
end
yasmin ismail
on 11 Jul 2024
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
on 13 Jul 2024
Edited: yasmin ismail
on 13 Jul 2024
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
on 13 Jul 2024
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
on 13 Jul 2024
Walter Roberson
on 13 Jul 2024
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
on 13 Jul 2024
x0 = [.7 .3 0 0
0 .5 .5 0
0 0 .8 .2
0 0 0 1]
But I suggest you first pass an introductory course to MATLAB free of costs before you continue working on your problem:
yasmin ismail
on 13 Jul 2024
So you want a random matrix with the above structure ?
d = rand(6,1);
du = 1- d;
x0 = diag([d;1]) + diag(du,1)
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 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)
% Reconstructed matrix
X = diag([x;1])+diag(1-x,1)
error = obj(x)
yasmin ismail
on 14 Jul 2024
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
on 14 Jul 2024
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)
% Reconstructed matrix
X = diag([x;1])+diag(1-x,1)
%Error
error = obj(x)
yasmin ismail
on 14 Jul 2024
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
on 15 Jul 2024
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.
p = [1 0 0 0 0 0 0];
R = [9; 8; 7; 6; 5; 4; 3];
CR = [0.45];
RHS1 = CR./(p*R)
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^(1/20)
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)
sol2 = nthroot(RHS1, 20)
Check by raising the solutions to the 20th power and comparing with RHS1.
format longg
[sol1^20; sol2^20; RHS1]
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 Comments
yasmin ismail
on 10 Jul 2024
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^2*initialState
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^2*initialState
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
on 10 Jul 2024
Edited: yasmin ismail
on 10 Jul 2024
Steven Lord
on 10 Jul 2024
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
on 11 Jul 2024
Categories
Find more on Mathematics in Help Center and File Exchange
on 10 Jul 2024
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
