5
$\begingroup$

I am trying to find an expression for the x[t+1] in the following recurrence relation

x[t + 1] == x[t]/P . x[t]

Where P is an invertible matrix with rational elements P[[i,j], and x[t] is a vector that starts as all ones.

Consider as a MWE:

n = 3;
P = HilbertMatrix[n];
x0 = ConstantArray[1, n];

RSolve[{x[t + 1] == x[t]/P . x[t], x[0] == x0}, 
 x[t] ∈ Vectors[n], t]

Is there any possibility to get an expression for x[t+1]?

An example with a little more applied background

I have a detection system with a binary circulant matrix A, a source vector x and a vector of observations y which are each independently Poisson-distributed with means equal to forward projection of the source through the matrix. As an example:

SeedRandom[123];
row1 = RandomChoice[{0, 1}, n];
A = NestList[RotateRight, row1, n - 1];

x = {10, 100, 10};

yMeans = A . x;

y = RandomVariate@*PoissonDistribution /@ yMeans;

Now define:

P = A\[Transpose] . A/A\[Transpose] . y;

Observe that P has the nice property when dotted with the least squares solution to A.xOLS ==y, it is a vector of all 1s:

xOLS = Inverse[A] . y;
P . xOLS

(*{1,1,1}*)

So the recurrence relation:

xGuess[t+1] == xGuess[t]/P.xGuess[t]

has a fixed point at the least-squares solution.

I don't want to go all the way to least-squares however because y is noisy, which is why I'm doing this iterative method, and I am interested in whether the recurrence relation xGuess[t+1] == xGuess[t]/P.xGuess[t] has a nice expression for xGuess[t+1].

Improve this question
$\endgroup$

2 Answers 2

Reset to default
5
$\begingroup$

With the parameters in the question and, for convenience,

Clear[x, y, z]
v[n_] := {x[n], y[n], z[n]}

both

RSolve[v[t + 1] == v[t]/P . v[t], v[t], t]

and

DSolveValue[D[v[t], t] == v[t]/P . v[t] - v[t], v[t], t]

return unevaluated. According to the documentation of RSolve, equations of the form v[t + 1] == v[t]/P . v[t] can be solved, if the rows of P are proportional to one another, but this clearly is not what the OP desires. So, let us turn to numerical solutions. Comparing

RecurrenceTable[{v[t + 1] == v[t]/P . v[t], v[0] == {1., 1., 1.}}, v[t], {t, 0, 40}];
pltr = ListLinePlot[Transpose[%], DataRange -> {0, 40}, 
    AxesLabel -> {t, "x,y,z"}, LabelStyle -> {10, Bold, Black}, 
    PlotLegends -> Placed[v[t], {.5, .5}]]

and its approximation,

NDSolveValue[{v'[t] == v[t]/P . v[t] - v[t], v[0] == {1, 1, 1}}, v[t], {t, 0, 10}];
pltd = Plot[%, {t, 0, 10}, PlotStyle -> Dashed]

we find,

Show[pltr, pltd]

enter image description here

where the solid curves are the recursion result and the dashed curves its ODE approximation. The two converge for t > 10. As A.Kato noted, the asymptotic solution has seven results.

Solve[v[t] == v[t]/P . v[t], v[t]];

However, at least one component of all but one of the solutions is zero. That one solution, which is what we seek, is given by

DeleteCases[%, a_ /; MemberQ[a, 0, Infinity]] // Flatten
(* {x[t] -> 7, y[t] -> 100, z[t] -> 14} *)

which agrees well with the curves.

Improve this answer
$\endgroup$
4
  • $\begingroup$ Thanks. FWIW, the nonzero stable point is also given by Inverse[P].ConstantArray[1,n] (since that makes the denominator all ones). I will keep trying to fiddle with the differential equation. AsymptoticDSolveValue can give a series expansion around t==0 to arbitrary degree, but I'm not really seeing a pattern yet in the series coefficients. $\endgroup$ Commented Jan 12 at 16:21
  • $\begingroup$ What level of accuracy are you seeking? $\endgroup$ Commented Jan 12 at 16:58
  • $\begingroup$ I'd like an exact expression for (ideally) the recurrence relation, and if not possible then for the differential equation. It's hard to explain here, but for my purposes I do not want the asymptotic fixed points, I want to stop the recurrence relation at a given condition that depends on the estimated noise in y. I currently just use NestWhile on the recurrence relation until the condition is met, but I wanted to see if I could derive an exact expression for the recurrence relation instead of having do it step-by-step. Sorry if this doesn't make a whole let of sense. $\endgroup$ Commented Jan 13 at 1:50
  • $\begingroup$ ... I know I'm leaving out a lot of details about my specific application, but I'm worried that would make the question more confusing. $\endgroup$ Commented Jan 13 at 1:50
4
$\begingroup$

The following does not answer your question, but hopefully might be of your help.

I don't know whether there exists a closed formula for x[t] or not. But you can compute first several terms quite easily:

ClearAll["Global`*"];
n = 3;
P = HilbertMatrix[n];
x0 = ConstantArray[1, n];
x[0] = x0;
x[t_Integer?Positive] := x[t] = x[t - 1]/P . x[t - 1];

Grid[Table[{k, x[k]}, {k, 0, 5}], Dividers -> All]

table

In order to find fixed points, let vec be a vector which represent a fixed point:

vec = Table[v[i], {i, n}]
(* {v[1], v[2], v[3]} *)

Then it has to satisfy:

vec == vec/P . vec

equ

Solve[%]

soln

Improve this answer
$\endgroup$
6
  • 1
    $\begingroup$ Thanks for your insights! I'm wondering if looking instead at the sort of analogous differential equation: $$ x'(t) = x(t) *\left(\frac{1}{P.x(t)} -1\right) $$ Would be helpful at all. soln[t_] = NDSolveValue[{v'[t] == v[t] (1/P . v[t] - 1), v[0] == x0}, v[t], {t, 0, 100}] But I doubt DSolve will be able to find an exact answer to this. $\endgroup$ Commented Jan 11 at 17:12
  • $\begingroup$ Please clarify the meaning of $\left(\frac{1}{P.x(t)}-1\right)$. The denominator $P.x(t)$ is a 3-component vector. $\endgroup$ Commented Jan 12 at 4:10
  • $\begingroup$ I mean element-wise division of each component of $P.x(t)$, if that makes sense $\endgroup$ Commented Jan 12 at 4:34
  • $\begingroup$ So if you define $\{z(t)_1,z(t)_2,z(t)_3\} = P.x(t)$ we have $\left(\frac{1}{P.x(t)} -1 \right) =\left\{ \frac{1}{z(t)_1} - 1,\frac{1}{z(t)_2} - 1, \frac{1}{z(t)_3}-1 \right\}$ $\endgroup$ Commented Jan 12 at 4:44
  • $\begingroup$ Well, if you want to change your problem from recursion to differential equation, I advise you to post it as another separate question :-) $\endgroup$ Commented Jan 12 at 7:50

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.