# Solve a system of M equations without specifying M

This post has been edited heavily in response to the comments: my initial question was very confusing.

I am trying to write my game-theoretic model in Wolfram Language. There is only one remaining step and I am not sure if it is within reach of Wolfram Language (or Mathematica). I describe a simplified version of my problem below.

# Cournot game

There are $M$ firms. Each firm $i$ can produce some output $q_i\geq0$ at a cost of $c_i$ per unit. Once the firms have (simultaneously) decided on their outputs, the price is set in the market according to an inverse demand function $P(Q)=\alpha-Q$, where $Q=q_1+\ldots+q_M$ is the aggregate output and $\alpha>0$ is market size parameter. Firm $i$'s profits are $q_iP(Q)-q_ic_i=q_i(\alpha-q_1-\ldots-q_M-c_i)$. Each firm maximizes its profits taking outputs of the opponents as given. By differentiating the profit function with respect to $q_i$ and imposing the necessary condition for maximizing a function, I obtain the equilibrium outputs: $$q_i^*=\frac{1}{2}\left(\alpha-\sum_{k\neq i}q_k^*-c_k\right) \text{ for } i=\{1,\ldots,M\}.$$ A quick derivation on the piece of paper gives me the explicit solution: $$q_i^*=\frac{1}{M+1}\left(\alpha-Mc_i+\sum_{k\neq i}c_k\right)$$ Here's the thing: I would like to obtain exactly this result (the above equation) in Mathematica. Is it possible? Note: I don't want to substitute any numbers for $M, c_1, \ldots, c_M$ or $\alpha$. They are general, exogenous parameters.

In other words, the above problem boils down to inverting a $M\times M$ matrix of the form: $$\left[\begin{array}{cccc} 2 & 1 & \ldots & 1\\ 1 & 2 & \ldots & 1\\ \vdots & \vdots & \ddots & \vdots \\ 1 & 1 & \ldots & 2 \end{array}\right].$$ Can one do this without specifying $M$?

I hope I expressed myself clearer this time.

-
If you have vectors c and q both of length m (any m) then Total[c q] does your calculation. You can create such vectors in many ways: ConstantArray[0,m] creates an m-length vector of zeros. RandomReal[{0,1},m] creates a vector of length m each of which is a random number in {0,1}. –  bill s Apr 22 at 17:48
No, that is not what I aim for. I have updated the description. I do not want to use specific values. –  Krzysztof Brzezinski Apr 22 at 18:00
I feel the closest you could get to a vector of symbolic length would be SymbolicTensors. However, they are mostly used in vector calculus like stuff (proving identities and so) and may not be what you are looking for. Perhaps it would help if you could elaborate a bit on what you want to do with them. –  Sjoerd C. de Vries Apr 22 at 18:11
I extended the description, I hope that clarifies things a bit. –  Krzysztof Brzezinski Apr 22 at 19:14
IdentityMatrix[M] - 1/(M + 1) ConstantArray[1, M]? –  kguler Apr 23 at 17:50

Maybe something like this:

cournotEq[m_] := (IdentityMatrix[m] -
1/(m + 1) ConstantArray[1, m]).Array[\[Alpha] - Subscript[c, #] &,  m]
cournotEq[2]//FullSimplify


$\left\{\frac{1}{3} \left(\alpha -2 c_1+c_2\right),\frac{1}{3} \left(\alpha +c_1-2 c_2\right)\right\}$

cournotEq[3]//FullSimplify


$\left\{\frac{1}{4} \left(\alpha -3 c_1+c_2+c_3\right),\frac{1}{4} \left(\alpha +c_1-3 c_2+c_3\right),\frac{1}{4} \left(\alpha +c_1+c_2-3 c_3\right)\right\}$

-
Almost. Can one make it work for cournotEq[M]? Without giving a concrete number? As it is, it is not possible. I want to get the exact equation I have given in the description (e.g. in a matrix form). –  Krzysztof Brzezinski Apr 23 at 20:59
@KrzysztofBrzezinski, I don't think it is possible to have it work for symbolic M. –  kguler Apr 23 at 21:23

You can invert the matrix $$\left[\begin{array}{cccc} 2 & 1 & \ldots & 1\\ 1 & 2 & \ldots & 1\\ \vdots & \vdots & \ddots & \vdots \\ 1 & 1 & \ldots & 2 \end{array}\right].$$ (for general m) using the Woodbury matrix identity. First, let one be the m-vector of all ones. Then the matrix you want to invert is expressible as

mat = IdentityMatrix[m] + Outer[Times, one, one];


To use the matrix inversion lemma (with the notation from the wikipedia article), let A = IdentityMatrix[m], U=one, V=Transpose[one], C=A; and the lemma says the inverse is:

matInv= IdentityMatrix[m] - Outer[Times, one, one]/(m + 1);


This holds for all m. We can check this for an specific m like so:

m = 5;
one = ConstantArray[1, m];
eye = IdentityMatrix[m];
mat = eye + Outer[Times, one, one];
matInv = eye - Outer[Times, one, one]/(m + 1);
mat.matInv == matInv.mat == eye


to which the response is True. To do this truly within Mathematica, one approach would be to incorporate knowledge of the inversion lemma.

-
I'm afraid this does not solve the problem. Can one make this method work without specifying m? –  Krzysztof Brzezinski Apr 24 at 17:10
If you want a numerical inverse to the matrix, then you must specify m. If you don't want a numerical inverse, then the formula for matInv is the answer: $I-1.1'/(m+1)$. –  bill s Apr 24 at 19:13