In God we trust, all others bring data.

– William Edwards Deming (1900-1993)

不知不觉 2020 年已经过去 1/12 了,1 月的文章也自然咕咕咕了(理直气壮.jpg)。

0. 一觉醒来变成了煤老板!

某天晚上,看着塞满的购物车,你嘀咕「好烦啊,下辈子想当家里有矿的帅哥!」。然而一觉醒来,看着身边堆积如山的煤炭,你不知道是快乐还是忧伤。

你看了看手上的劳力士,记忆如潮水般袭来。你叫朱一蛋,是西虹市最大的煤矿公司老总,正准备回公司开晨会。哎,有钱人的快乐,往往就是这么简简单单……

秘书突然出现,并打断了你的回忆:“老板,晨会时间快到了。”你点了点头,和秘书回到了公司会议室。拿起桌上的财务报表,你大致了解了公司的情况:你的煤炭公司和西虹市的电力公司、铁路运输公司是当地经济三巨头,而且三家公司紧密相连。

  • 煤炭公司每生产价值 1 元的煤炭,需要消耗价值 0.25 元的电和价值 0.35 元的铁路运输服务;
  • 电力公司每生产价值 1 元的电,需要消耗价值 0.40 元的煤炭、价值 0.05 元的电和价值 0.10 元的运输服务;
  • 铁路运输公司每提供价值 1 元的运输服务,需要消耗价值 0.45 元的煤炭、价值 0.10 元的电和价值 0.10 元的运输服务。

晨会开始了,销售部长发言,表示上个月拿到了价值 10000 元的煤炭订单,需要在这个月交付完成。生产部长发言,请求你下达生产数量指示。

心不在焉的你随口回答:“那就生产 10000 元煤炭吧。”

话音刚落,你看到众人怀疑的神色,小声问秘书:“一个月生产 10000 元的煤炭有问题吗?”

秘书:“没有,但是……”

你:“但是什么?”

秘书:“但是生产 10000 元的煤炭,需要消耗 2500 元的电和 3500 元的铁路运输服务。”

你:“没问题啊,花钱买就行了。”

秘书:“而电力公司生产 2500 元的电,需要消耗 1000 元的煤炭、125 元的电和 250 元的运输服务。”

你:“哦,那就再加 1000,生产 11000 元就够了吧。”

秘书:“运输公司提供 3500 元的运输服务,需要消耗 1575 元的煤炭、350 元的电和 350 元的运输服务,加起来是 2575 元的煤炭、475 元的电和 600 元的运输服务。”

你:“所以一共是 12575 元的煤炭?”

秘书:“而这 475 元的电,需要消耗 190 元的煤炭、23.75 元的电和 47.5 元的运输服务。”

……

5 分钟过去了,然而秘书的计算完全没有停下来的迹象,你意识到事情似乎并不简单!

1. 消耗系数方阵(Input Coefficient Matrix)

我们不妨用矩阵的形式表示上面的计算过程:

\left\{ \begin{array}{l}
  a_0 = \left(\begin{array}{c}
    10000\\
    0\\
    0
  \end{array}\right)\\
  a_1 = \left(\begin{array}{ccc}
    0 & 0.4 & 0.45\\
    0.25 & 0.05 & 0.1\\
    0.35 & 0.1 & 0.1
  \end{array}\right) \left(\begin{array}{c}
    10000\\
    0\\
    0
  \end{array}\right) = \left(\begin{array}{c}
    0\\
    2500\\
    3500
  \end{array}\right)\\
  a_2 = \left(\begin{array}{ccc}
    0 & 0.4 & 0.45\\
    0.25 & 0.05 & 0.1\\
    0.35 & 0.1 & 0.1
  \end{array}\right) \left(\begin{array}{c}
    0\\
    2500\\
    3500
  \end{array}\right) = \left(\begin{array}{c}
    2575\\
    475\\
    600
  \end{array}\right)\\
  \ldots\\
  a_n = \left(\begin{array}{ccc}
    0 & 0.4 & 0.45\\
    0.25 & 0.05 & 0.1\\
    0.35 & 0.1 & 0.1
  \end{array}\right) a_{n - 1} = \left(\begin{array}{ccc}
    0 & 0.4 & 0.45\\
    0.25 & 0.05 & 0.1\\
    0.35 & 0.1 & 0.1
  \end{array}\right)^n a_0\\
  A = \left(\begin{array}{ccc}
    0 & 0.4 & 0.45\\
    0.25 & 0.05 & 0.1\\
    0.35 & 0.1 & 0.1
  \end{array}\right) \Rightarrow a_n = A^n a_0\\
  a = \lim\limits_{n \rightarrow + \infty} a_0 + a_1 + \cdots + a_n = \lim\limits_{n
  \rightarrow + \infty} (I + A + A^2 + \ldots + A^n) a_0
\end{array} \right.

可以留意到对任意给定的产出 x,Ax 就是对应的投入(消耗),因此矩阵 A 被称为 消耗系数方阵(Input Coefficient Matrix)。那么 a 到底是多少呢?回忆等比数列的求和公式,易得:

(1 - q) (1 + q + q^2 + \cdots + q^n) = 1 - q^n

容易推广到矩阵的情况:

(I - A) (I + A + A^2 + \cdots + A^n) = I - A^n

回忆等比数列的性质,如果 q 的 n 次方趋近于 0,那么和趋近 1/(1-q):

\begin{array}{rl}
  & \lim\limits_{n \rightarrow + \infty} A^n = 0\\
  \Rightarrow \; a = & \lim\limits_{n \rightarrow + \infty} (I + A + A^2 + \ldots +
  A^n) a_0\\
  = & \lim\limits_{n \rightarrow + \infty} (I - A)^{- 1} (I - A^n) a_0\\
  = & (I - A)^{- 1} \lim\limits_{n \rightarrow + \infty} (I - A^n) a_0\\
  = & (I - A)^{- 1} a_0
\end{array}

为了简化讨论,我们引入如下的定义:

  • 如果矩阵 A 每个元素都大于等于 0,则称矩阵 A 非负(Nonnegative),向量同理;
  • 如果对任意的非负向量 a0,解出的向量 a 都是非负的,则称矩阵 A 可行。(任何非负外部需求有相应非负的总产值,称为合理的模型。)

什么样的非负矩阵 A 是可行的呢?显然首先需要 A^n->0,以下的矩阵均满足这一条件1

  • A 的行和最大值 \max \left( \sum_{k = 1}^n A_{1 k}, \sum_{k = 1}^n A_{2 k}, \ldots, \sum_{k = 1}^n A_{\tmop{nk}} \right) 小于 1;
  • A 的列和最大值 \max \left( \sum_{k = 1}^n A_{k 1}, \sum_{k = 1}^n A_{k 2}, \ldots, \sum_{k = 1}^n A_{\tmop{kn}} \right) 小于 1;
  • A 的所有特征值绝对值的最大值 \max (| \lambda_1 |, | \lambda_2 |, \ldots, | \lambda_n |) 小于 1(充要条件)2

一般对于消耗系数方阵 A 而言,列和最大值都会小于 1,否则该列对应的商品每生产一件都在亏钱(毕竟这不是péi总的小说😂),所以这个条件一般都能满足。

接下来我们要使 a 可解,显然只要 I-A 可逆即可满足。

最后对任意的非负向量 a0,解出来的 a 也是非负向量。显然如果矩阵 (I-A)^-1 非负,即可保证成立。3

解决了可解性的问题,现在我们回头再看,(I-A)^-1 有什么实际的意义?回顾消耗系数方阵的定义,对每一个产出 X 而言,AX 都是其投入,于是增量就是产出减去投入:Y = X-AX =(I-A)X。我们实际求得的,是在保证生产系统稳定运行下,满足增量需求的最低产量。

想明白这一切的你,高兴地打开 Octave(MATLAB/Excel/R/etc):

A = [0 0.4 0.45;0.25 0.05 0.1;0.35 0.1 0.1];
Y = [10000;0;0];
X = (eye(3) - A)\Y
# X =

   14565.82633
    4481.79272
    6162.46499

顺畅安排好生产任务后,你度过了无忧无虑的第一个月。哎,有钱人的快乐,往往就是这么简简单单……

2. 又好又快,均衡发展

顺利完成第一个月生产任务的你,收到了市政府的会议文件。文件指出,要坚持以经济建设为中心的基本原则,同时也要保证经济的均衡发展。相似的情节再次展开,生产部长又来请求生产任务了……

你揉了揉脑袋,「以经济建设为中心」暂时想不出来,那么怎么才能保证经济的均衡发展呢?你很快意识到,首先需要衡量增长的快慢,然后让三家公司的增长率完全一致就可以了!

由于已经知道产出向量是 X,投入向量是 AX,增长率自然就是 X 的每一个元素,除以 AX 的对应元素。「咦,怎么有点熟悉的感觉?唉,如果 A 不是矩阵而是个常数就好了……等等,把矩阵当作常数?!」

灵光一闪,你喊出了三个字:“特!征!值!”

如果矩阵 A 恰好有一个小于 1 的正特征值 l,而且这个特征值对应的特征向量也是正的,只要令 X 为特征向量,则 AX=lX,增长率自然都是 1/l。于是你又美滋滋打开了 Octave:

[V, L] = eig(A)
# V =

   0.704192   0.814191  -0.031341
   0.420672  -0.326962  -0.746050
   0.571969  -0.479781   0.665152

# L =

   0.604458          0          0
          0  -0.425804          0
          0          0  -0.028654

简单验算后你发现真的存在:

\begin{array}{l}
  Av_1 = \left(\begin{array}{ccc}
    0 & 0.4 & 0.45\\
    0.25 & 0.05 & 0.1\\
    0.35 & 0.1 & 0.1
  \end{array}\right) \left(\begin{array}{c}
    0.704192\\
    0.420672\\
    0.571969
  \end{array}\right) = \left(\begin{array}{c}
    0.42565\\
    0.25428\\
    0.34573
  \end{array}\right)\\
  \lambda_1 v_1 = 0.604458 \cdot \left(\begin{array}{c}
    0.704192\\
    0.420672\\
    0.571969
  \end{array}\right) = \left(\begin{array}{c}
    0.42565\\
    0.25428\\
    0.34573
  \end{array}\right) = Av_1
\end{array}

于是只要让这一个月的投入尽可能接近特征向量的倍数即可。上个月的结余为 (4566, 4482, 6162),故可以令投入 AX = 6484 * v1 = (4566, 2728, 3709),因此产出 X = 1/l * (4566, 2728, 3709) = (7553, 4513, 6136),投入产出比 1/l = 1.6544。

完成与电力公司和运输公司的生产协商后,你抬了抬手上的劳力士,又是枯燥的一个月呢……

3. 多快好省真的存在吗?

虽然对矩阵 A 取得了阶段性的胜利,但这个方法是可行的吗?具体而言:

  1. 什么样的矩阵才能保证拥有正的特征值和特征向量?
  2. 这样求出的投入产出比是最大的吗(经济增长最快)?

对于问题 1,我们有著名的 Perron-Frobenius 定理:一个不可约的非负方阵,特征值绝对值的最大值一定为正,而且它对应的(左/右)特征向量只有一个,且这个特征向量是正的。虽然大部分的矩阵都是不可约的,但是可约的矩阵也不代表不存在正的特征值与特征向量,只是也许存在且好几个也许不存在罢了。

写这篇文章之前纠结了好久要不要把这个定理的证明写出来,在翻阅了三个不同的证明后放弃了这个想法……虽然这个定理很美妙,但是证明的完整篇幅比这篇文章的几倍还长😂简要叙述一下其中一个证明的思路4:类似相似标准型的做法(但还是很巧妙!),不可约的矩阵最后都会置换相似于一个行和相等的矩阵(称为标准型),假设行和为 q,很容易证明 q 就是矩阵的特征值,同时置换矩阵(全是正数)给出了一个特征向量。我就知道肯定没人会看的233

而问题 2 取决于衡量快慢的标准,如果只需要让增长率的均值(总和)最大化,像开头一样只生产一项就会「更快」,但显然是不合理的,没有消费的情况下不可能持续下去。可以证明「正特征向量法」能使增长率的最小值最大化(让木桶的短板最长),证明也跳过好了……

4. 哎呀,崩溃啦!

有了「正特征向量法」后,你已经几个月没去过公司了,因为如果 X 是 A 的特征向量,则 A^2 X = A (AX) = A (\lambda X) = \lambda (AX) = \lambda^2 X ,也就是 X 也是 A^2, A^3, … A^n 乃至 A^-1, A^-2, …, A^-n 的特征向量,在这几个月里三巨头都在蓬勃发展,而且似乎永不停歇……

然而今天,秘书突然发来消息,铁路运输公司快要倒闭了!

你心想,完全不可能!正特征向量法保证了每次增长率都会是正数 1/l,怎么可能就变成负值了呢。

然而打开 Octave 一算,竟然是真的:

X = [4566;2728;3709]
A = [0 0.4 0.45;0.25 0.05 0.1;0.35 0.1 0.1];
B = A^-1;
(B^4)*X
# ans =

    38461.91922
   121962.31836
   -62727.86294

“怎么会这样?这根本不科学!”

几天过后,煤炭公司也倒闭了,监狱里也从此多了个疯子,每天都在自言自语:「为什么会变成这样……」

一周目结束,恭喜你获得了 Bad End😂


十分让人意外的是,对于偏离正特征向量的投入 X(哪怕只是因为计算误差偏了一点点),总会存在一个 n,使得 B^n * X 中至少有一项小于 0!5 这意味着对经济发展而言,不平衡才是常态,平衡的概率,仅仅是 n 维空间里的一条射线,无限接近于 0。6

如何才能延长崩溃的时间呢?以下结论不难看出:

  • 定期调整投入向量,使其尽可能接近正特征向量;
  • 优化产业结构,改变 A 的正特征向量方向;
  • 促进消费,显然如果把每次的产出增量都消耗完,经济永远不会崩溃;
  • ……

虽然出人意料,但仔细想想又在情理之中,以陈木法先生的评论作为结语吧:

上述结论甚妙,它表明:依据所具备的客观条件,经济有它自身的合理的发展速度。太快了不好,太慢了也不行。其实点破了也就很容易明白。我们可以在一夜之间办起无数的工厂,但没有足够的电力和交通,这些厂子能运转起来吗?我们可创办很多很多的大学,但没有教师,没有足够财力的支持,资料室是空的,实验仪器一样也没有,能算得上大学吗?7

5. 今回はここまで

到这里就结束啦~
到这里就结束啦~

5.1 历史

投入产出法的诞生有一个有趣的小插曲。最早提出投入产出法的是哈佛大学教授 Leontief,当时他和同事希望根据各种统计表计算出美国各行业的生产总值,然而这需要计算一个 24 阶的逆矩阵(但愿你能回想起计算逆矩阵的痛苦)。面对这一矩阵他们犯了难,根据他们估计,即使是一周工作 7 天,算出结果也要花上几百年的时间。恰好哈佛大学发明了第一台原始计算机,于是他们租借这台计算机完成了这一本不可能完成的计算,Leontief 也因此获得了 1973 年的诺贝尔经济学奖。

然而事后他们希望为此付费的时候却被劳工部门阻止了,原因是政府部门有一项政策,货物可以购买,而服务不能购买。后来他们让劳工统计局向哈佛大学开出一张购买固定资产的订单,绕过了这一限制,而这项固定资产叫“逆矩阵”。8

文中提及的正特征向量法和崩溃定理,都是由我国著名数学家华罗庚先生提出的9。尽管最早是用来证明计划经济的可行性,这一方法很快就推广到市场经济的情形。在随机分析方面的处理,陈木法先生也有贡献7

5.2 现状与局限性

投入产出法的局限性是很明显的:

  • 没有考虑生产能力的限制,对任意的正需求都有正产出;
  • 生产能力不是可持续的,比如煤矿越往后越难挖;
  • 生产能力会被偶然性因素影响,比如农业收成严重依赖天气;
  • 模型必然有误差,然而误差是致命的,会加速经济的崩溃;
  • ……

然而尽管如此,作为一种最简单实用的模型,投入产出法早在 1968 年,就被联合国列为国民经济核算的工具。我国从 1974 年开始编制投入产出表,每隔 2,3 年统计一次,数据均可在国家统计局网站下载10。感兴趣的也可以试试用正特征向量法分析~

啊,差点忘了,投入产出表要转化成消耗系数方阵,只要(中间产物的部分)每一列除以列和(行和)就行了~(为什么?)

呼,去睡了,好困……

  1. 也许你已经留意到,要定义矩阵的极限,需要先定义矩阵的范数,这三个条件实际上是矩阵的三种由欧氏范数诱导的 算子范数(Operator Norm),满足 norm(AB) <= norm(A)norm(B),故若 norm(A) < 1,则 norm(A^n) <= norm(A)^n -> 0。 

  2. 也被称为 谱半径(Spectral Radius)。必要性是显然的,如果 lim A^n = 0, x 为任一特征向量,l 为对应的特征值,则 Ax = \lambda x \Rightarrow \lim_{k \rightarrow + \infty} \lambda^k X = (\lim_{k \rightarrow + \infty} A^k) X = 0 \Rightarrow | \lambda | < 1 ,充分性的证明需要使用 Jordan 标准型。 

  3. 矩阵 (I-A)^-1 被称为 Leontief 逆矩阵。 

  4. 华罗庚《高等数学引论余篇》。 

  5. 学过概率论的可能会联想到马尔可夫链的极限定理,两者确实是有联系的~ 

  6. 按测度论的观点就是 0。 

  7. 《随机过程导论》及后面的参考文献。  2

  8. 《女士品茶(The Lady Tasting Tea)》,Chapter 17。 

  9. 《计划经济大范围最优化的数学理论》,I-X。 

  10. 最近的一次是 2019 年发布的 2017 年投入产出表: http://data.stats.gov.cn/ifnormal.htm?u=/files/html/quickSearch/trcc/trcc01.html&h=740 。