ISTA与稀疏编码

本人习惯

注释

我自己的想法或者例子，不必要

代码

function [...] = fun(...)
....
end

稀疏编码

稀疏编码是什么鬼

图像可以认为是n×m×3n \times m \times 3n×m×3的一个矩阵，随着图像清晰度的提升，以目前主流的1920*1800图像为例，单单一个图片就可能有1920×1080×3≥6×1061920 \times 1080 \times 3 \geq 6 \times 10^{6}1920×1080×3≥6×106,在假设每个位置使用256色(282^828种颜色需要8bit，也就是1B)，一张图片的大小已经大于6MB，未经压缩的1080P视频每分钟数据量约为1920×1080×3×30×60≥1×10101920 \times 1080 \times 3 \times 30 \times 60 \geq 1 \times 10^{10}1920×1080×3×30×60≥1×1010,一分钟超过一个G的存储量显然目前的存储技术无法解决，只能从无损压缩图像->无损压缩视频这个思路出发。虽然JPEG提供了一种无损压缩算法，但是速度较慢。一般情况下，对于实时性要求高或数据传输成本高的情况(例如医疗影像\卫星图像)，有损压缩成为了一种必要之恶。
稀疏编码是一种依托优化思想出现的压缩方式，他从两个方向出发对图像进行压缩：

保真度
对于原始的信号X，我们希望有一个信号Y充分逼近X，有如下关系：
X=Y+ϵ\begin{aligned} X = Y + \epsilon \end{aligned} X=Y+ϵ
其中ϵ\epsilonϵ为噪声，我们在这里希望噪声造成的影响最小，即：
Y∗=arg⁡min⁡YB(X,Y)\begin{aligned} Y^* = {\arg\min}_Y B(X,Y) \end{aligned} Y∗=argminYB(X,Y)
其中B(X,Y)=∣∣Y−X∣∣pp=∣∣ϵ∣∣ppB(X,Y)=||Y-X||_p^p = ||\epsilon||_p^pB(X,Y)=∣∣Y−X∣∣pp=∣∣ϵ∣∣pp
从B的定义考虑，若我们选择p=2，B(X,Y)可以认为是均方误差（MSE），我们希望其误差尽量小。也是一种合理的解释呢！
接下来，我们给定一个Y的表示方法
Y=D⋅Z\begin{aligned} Y = D \cdot Z \end{aligned} Y=D⋅Z
显然，如果D是n维空间的一组标准正交基，Z就相当于是Y在D上映射的坐标，且对于指定的D，Z有唯一表示，我们希望Z有某种优秀的性质。那么我们有几个方向可以考虑
- 推断D的结构，通过数据对D进行优化，保证得到的Z拥有好的性质，但是这太难了，对于不同的数据，很难用矩阵论的方法对D进行推断。
- 如果上述结构不进行推断，仍可以通过机器学习算法例如BP神经网络优化来优化这个矩阵D，使得D对于一个特定场景具有较好的性质。确实是一种好的思想。但是如果D发生变化，Z将无法通过简单的矩阵乘法变回信号Y，显然如果空间的基发送变化，这组基上的坐标都将失去原来的意义。需要大量的数据才能保证在某一问题上D有良好的效果。只能考虑继续优化思路。
- 我们先给定一组标准正交基S，再给定一些列向量，写成矩阵P，如果D = [S P], 那么P中所有向量均可由S线性表出，此时D(之后称为字典)可以称为一个超完备的集合，此时Z的表示不再唯一，为后续对于Z某种性质的优化提供了可能性。(*对于ISTA我们考虑到这里就结束了，但是仍然可以对此进行改进，如果有好的想法欢迎私信探讨)
```
关于字典，我觉得是他通过一个拥有优秀性质的Z表示了蕴含复杂信息的Y
```
稀疏性
对于一个向量Z，Z中非0元素个数为u，则称Z为u稀疏的(u-Sparse)。对于一个充分稀疏的向量Z，使用Z{Site, Value}表示更加节省空间(类似图论算法中的两种存图结构：邻接矩阵和邻接表)，是一种可行的压缩方式。
```
Z = [0,0,1,0,0,0,3,0,0,0,0,0,1,0,0,2,0,5,0,0,0,0,0,0,0,0,0,0,0]
Z{Site, Value} = [{2,1},{6,3},{12,1},{15,2},{17,5}]
很明显sizeof(Z)>sizeof(Z{Site, Value})
```
那么如何对通过优化保证稀疏性呢？
Z∗=arg⁡min⁡ZS(X,Y)=∣∣Z∣∣0\begin{aligned} Z^* = {\arg\min}_Z S(X,Y) = ||Z||_0 \end{aligned} Z∗=argminZS(X,Y)=∣∣Z∣∣0
是一种很好的想法，但是L-0范数显然是非凸的，优化难度巨大，属于NP-Hard问题，只能智取，不可强攻。
IDEAs:
1. 使用贪心算法(one can say it nearly always choose all the suboptimal solution)
2. 使用群体智能算法(some said only works well on low-dimension question)
3. 放宽条件(some called of ’松弛‘， which may be a good idea)
根据大佬的证明，将0范数放宽到1范数是可行的，故原问题转化为：
Z∗=arg⁡min⁡ZS(Z)=∣∣Z∣∣1\begin{aligned} Z^* = {\arg\min}_Z S(Z) = ||Z||_1 \end{aligned} Z∗=argminZS(Z)=∣∣Z∣∣1

优秀的糅合怪LASSO格式的稀疏优化

我们有两种思路将上述问题糅合起来，就是稀疏编码的目标：
在有足够保真度的区域找到稀疏性最好的点：
Z∗=arg⁡min⁡ZS(Z)\begin{aligned} Z^* = {\arg\min}_Z S(Z) \end{aligned} Z∗=argminZS(Z) s.t.B(X,D⋅Z)<=e\begin{aligned} s.t. B(X,D \cdot Z) <= e \end{aligned} s.t.B(X,D⋅Z)<=e
在稀疏性足够好的区域找到保真度最高的点：Z∗=arg⁡min⁡ZB(X,D⋅Z)\begin{aligned} Z^* = {\arg\min}_Z B(X,D \cdot Z) \end{aligned} Z∗=argminZB(X,D⋅Z) s.t.S(Z)<=u\begin{aligned} s.t.S(Z) <= u \end{aligned} s.t.S(Z)<=u
使用罚函数法把上述两个问题都可以转化为如下LASSO格式：Z∗=arg⁡min⁡ZB(X,D⋅Z)+λS(Z)\begin{aligned} Z^* = {\arg\min}_Z B(X,D \cdot Z) + \lambda S(Z) \end{aligned} Z∗=argminZB(X,D⋅Z)+λS(Z)
这就是我们接下来要优化的问题了。

迭代软阈值算法

软阈值函数

hθ(x)=sign(x)×max(abs(x)−θ,0)h_\theta(x) = sign(x) \times max(abs(x)-\theta, 0)hθ(x)=sign(x)×max(abs(x)−θ,0)

Other`s Description about Soft Thresholding： https://blog.csdn.net/jbb0523/article/details/52103257

ISTA

由LASSO问题的解形式可以得到(Z0=0⃗Z_0=\vec 0Z0=0)：
Zt+1=Zt−1LWdT(WdZt−X)=1LWdTX+(I−1LWdTWd)Zt\begin{aligned} Z_{t+1} = Z_t - \frac{1}{L}W_d^T(W_dZ_t-X) \\ = \frac{1}{L}W_d^TX+(I-\frac{1}{L}W_d^TW_d)Z_t \end{aligned} Zt+1=Zt−L1WdT(WdZt−X)=L1WdTX+(I−L1WdTWd)Zt
若取
S=I−1LWdTWdWe=1LWdTL≥rank(WdTWd)S = I - \frac{1}{L}W_d^TW_d\\ W_e = \frac{1}{L}W_d^T\\ L \ge rank(W_d^TW_d) S=I−L1WdTWdWe=L1WdTL≥rank(WdTWd)
则迭代公式可化为：
Zt+1=WeX+SZt\begin{aligned} Z_{t+1} = W_eX+SZ_t \end{aligned} Zt+1=WeX+SZt
若每一次迭代后经过一次激活则有：
Zt+1=hθ(WeX+SZt)\begin{aligned} Z_{t+1} = h_\theta(W_eX+SZ_t) \end{aligned} Zt+1=hθ(WeX+SZt)
则有如下算法模型（RNN）：

#mermaid-svg-4vuN1MjXmBoSUadg .label{font-family:'trebuchet ms', verdana, arial;font-family:var(--mermaid-font-family);fill:#333;color:#333}#mermaid-svg-4vuN1MjXmBoSUadg .label text{fill:#333}#mermaid-svg-4vuN1MjXmBoSUadg .node rect,#mermaid-svg-4vuN1MjXmBoSUadg .node circle,#mermaid-svg-4vuN1MjXmBoSUadg .node ellipse,#mermaid-svg-4vuN1MjXmBoSUadg .node polygon,#mermaid-svg-4vuN1MjXmBoSUadg .node path{fill:#ECECFF;stroke:#9370db;stroke-width:1px}#mermaid-svg-4vuN1MjXmBoSUadg .node .label{text-align:center;fill:#333}#mermaid-svg-4vuN1MjXmBoSUadg .node.clickable{cursor:pointer}#mermaid-svg-4vuN1MjXmBoSUadg .arrowheadPath{fill:#333}#mermaid-svg-4vuN1MjXmBoSUadg .edgePath .path{stroke:#333;stroke-width:1.5px}#mermaid-svg-4vuN1MjXmBoSUadg .flowchart-link{stroke:#333;fill:none}#mermaid-svg-4vuN1MjXmBoSUadg .edgeLabel{background-color:#e8e8e8;text-align:center}#mermaid-svg-4vuN1MjXmBoSUadg .edgeLabel rect{opacity:0.9}#mermaid-svg-4vuN1MjXmBoSUadg .edgeLabel span{color:#333}#mermaid-svg-4vuN1MjXmBoSUadg .cluster rect{fill:#ffffde;stroke:#aa3;stroke-width:1px}#mermaid-svg-4vuN1MjXmBoSUadg .cluster text{fill:#333}#mermaid-svg-4vuN1MjXmBoSUadg div.mermaidTooltip{position:absolute;text-align:center;max-width:200px;padding:2px;font-family:'trebuchet ms', verdana, arial;font-family:var(--mermaid-font-family);font-size:12px;background:#ffffde;border:1px solid #aa3;border-radius:2px;pointer-events:none;z-index:100}#mermaid-svg-4vuN1MjXmBoSUadg .actor{stroke:#ccf;fill:#ECECFF}#mermaid-svg-4vuN1MjXmBoSUadg text.actor>tspan{fill:#000;stroke:none}#mermaid-svg-4vuN1MjXmBoSUadg .actor-line{stroke:grey}#mermaid-svg-4vuN1MjXmBoSUadg .messageLine0{stroke-width:1.5;stroke-dasharray:none;stroke:#333}#mermaid-svg-4vuN1MjXmBoSUadg .messageLine1{stroke-width:1.5;stroke-dasharray:2, 2;stroke:#333}#mermaid-svg-4vuN1MjXmBoSUadg #arrowhead path{fill:#333;stroke:#333}#mermaid-svg-4vuN1MjXmBoSUadg .sequenceNumber{fill:#fff}#mermaid-svg-4vuN1MjXmBoSUadg #sequencenumber{fill:#333}#mermaid-svg-4vuN1MjXmBoSUadg #crosshead path{fill:#333;stroke:#333}#mermaid-svg-4vuN1MjXmBoSUadg .messageText{fill:#333;stroke:#333}#mermaid-svg-4vuN1MjXmBoSUadg .labelBox{stroke:#ccf;fill:#ECECFF}#mermaid-svg-4vuN1MjXmBoSUadg .labelText,#mermaid-svg-4vuN1MjXmBoSUadg .labelText>tspan{fill:#000;stroke:none}#mermaid-svg-4vuN1MjXmBoSUadg .loopText,#mermaid-svg-4vuN1MjXmBoSUadg .loopText>tspan{fill:#000;stroke:none}#mermaid-svg-4vuN1MjXmBoSUadg .loopLine{stroke-width:2px;stroke-dasharray:2, 2;stroke:#ccf;fill:#ccf}#mermaid-svg-4vuN1MjXmBoSUadg .note{stroke:#aa3;fill:#fff5ad}#mermaid-svg-4vuN1MjXmBoSUadg .noteText,#mermaid-svg-4vuN1MjXmBoSUadg .noteText>tspan{fill:#000;stroke:none}#mermaid-svg-4vuN1MjXmBoSUadg .activation0{fill:#f4f4f4;stroke:#666}#mermaid-svg-4vuN1MjXmBoSUadg .activation1{fill:#f4f4f4;stroke:#666}#mermaid-svg-4vuN1MjXmBoSUadg .activation2{fill:#f4f4f4;stroke:#666}#mermaid-svg-4vuN1MjXmBoSUadg .mermaid-main-font{font-family:"trebuchet ms", verdana, arial;font-family:var(--mermaid-font-family)}#mermaid-svg-4vuN1MjXmBoSUadg .section{stroke:none;opacity:0.2}#mermaid-svg-4vuN1MjXmBoSUadg .section0{fill:rgba(102,102,255,0.49)}#mermaid-svg-4vuN1MjXmBoSUadg .section2{fill:#fff400}#mermaid-svg-4vuN1MjXmBoSUadg .section1,#mermaid-svg-4vuN1MjXmBoSUadg .section3{fill:#fff;opacity:0.2}#mermaid-svg-4vuN1MjXmBoSUadg .sectionTitle0{fill:#333}#mermaid-svg-4vuN1MjXmBoSUadg .sectionTitle1{fill:#333}#mermaid-svg-4vuN1MjXmBoSUadg .sectionTitle2{fill:#333}#mermaid-svg-4vuN1MjXmBoSUadg .sectionTitle3{fill:#333}#mermaid-svg-4vuN1MjXmBoSUadg .sectionTitle{text-anchor:start;font-size:11px;text-height:14px;font-family:'trebuchet ms', verdana, arial;font-family:var(--mermaid-font-family)}#mermaid-svg-4vuN1MjXmBoSUadg .grid .tick{stroke:#d3d3d3;opacity:0.8;shape-rendering:crispEdges}#mermaid-svg-4vuN1MjXmBoSUadg .grid .tick text{font-family:'trebuchet ms', verdana, arial;font-family:var(--mermaid-font-family)}#mermaid-svg-4vuN1MjXmBoSUadg .grid path{stroke-width:0}#mermaid-svg-4vuN1MjXmBoSUadg .today{fill:none;stroke:red;stroke-width:2px}#mermaid-svg-4vuN1MjXmBoSUadg .task{stroke-width:2}#mermaid-svg-4vuN1MjXmBoSUadg .taskText{text-anchor:middle;font-family:'trebuchet ms', verdana, arial;font-family:var(--mermaid-font-family)}#mermaid-svg-4vuN1MjXmBoSUadg .taskText:not([font-size]){font-size:11px}#mermaid-svg-4vuN1MjXmBoSUadg .taskTextOutsideRight{fill:#000;text-anchor:start;font-size:11px;font-family:'trebuchet ms', verdana, arial;font-family:var(--mermaid-font-family)}#mermaid-svg-4vuN1MjXmBoSUadg .taskTextOutsideLeft{fill:#000;text-anchor:end;font-size:11px}#mermaid-svg-4vuN1MjXmBoSUadg .task.clickable{cursor:pointer}#mermaid-svg-4vuN1MjXmBoSUadg .taskText.clickable{cursor:pointer;fill:#003163 !important;font-weight:bold}#mermaid-svg-4vuN1MjXmBoSUadg .taskTextOutsideLeft.clickable{cursor:pointer;fill:#003163 !important;font-weight:bold}#mermaid-svg-4vuN1MjXmBoSUadg .taskTextOutsideRight.clickable{cursor:pointer;fill:#003163 !important;font-weight:bold}#mermaid-svg-4vuN1MjXmBoSUadg .taskText0,#mermaid-svg-4vuN1MjXmBoSUadg .taskText1,#mermaid-svg-4vuN1MjXmBoSUadg .taskText2,#mermaid-svg-4vuN1MjXmBoSUadg .taskText3{fill:#fff}#mermaid-svg-4vuN1MjXmBoSUadg .task0,#mermaid-svg-4vuN1MjXmBoSUadg .task1,#mermaid-svg-4vuN1MjXmBoSUadg .task2,#mermaid-svg-4vuN1MjXmBoSUadg .task3{fill:#8a90dd;stroke:#534fbc}#mermaid-svg-4vuN1MjXmBoSUadg .taskTextOutside0,#mermaid-svg-4vuN1MjXmBoSUadg .taskTextOutside2{fill:#000}#mermaid-svg-4vuN1MjXmBoSUadg .taskTextOutside1,#mermaid-svg-4vuN1MjXmBoSUadg .taskTextOutside3{fill:#000}#mermaid-svg-4vuN1MjXmBoSUadg .active0,#mermaid-svg-4vuN1MjXmBoSUadg .active1,#mermaid-svg-4vuN1MjXmBoSUadg .active2,#mermaid-svg-4vuN1MjXmBoSUadg .active3{fill:#bfc7ff;stroke:#534fbc}#mermaid-svg-4vuN1MjXmBoSUadg .activeText0,#mermaid-svg-4vuN1MjXmBoSUadg .activeText1,#mermaid-svg-4vuN1MjXmBoSUadg .activeText2,#mermaid-svg-4vuN1MjXmBoSUadg .activeText3{fill:#000 !important}#mermaid-svg-4vuN1MjXmBoSUadg .done0,#mermaid-svg-4vuN1MjXmBoSUadg .done1,#mermaid-svg-4vuN1MjXmBoSUadg .done2,#mermaid-svg-4vuN1MjXmBoSUadg .done3{stroke:grey;fill:#d3d3d3;stroke-width:2}#mermaid-svg-4vuN1MjXmBoSUadg .doneText0,#mermaid-svg-4vuN1MjXmBoSUadg .doneText1,#mermaid-svg-4vuN1MjXmBoSUadg .doneText2,#mermaid-svg-4vuN1MjXmBoSUadg .doneText3{fill:#000 !important}#mermaid-svg-4vuN1MjXmBoSUadg .crit0,#mermaid-svg-4vuN1MjXmBoSUadg .crit1,#mermaid-svg-4vuN1MjXmBoSUadg .crit2,#mermaid-svg-4vuN1MjXmBoSUadg .crit3{stroke:#f88;fill:red;stroke-width:2}#mermaid-svg-4vuN1MjXmBoSUadg .activeCrit0,#mermaid-svg-4vuN1MjXmBoSUadg .activeCrit1,#mermaid-svg-4vuN1MjXmBoSUadg .activeCrit2,#mermaid-svg-4vuN1MjXmBoSUadg .activeCrit3{stroke:#f88;fill:#bfc7ff;stroke-width:2}#mermaid-svg-4vuN1MjXmBoSUadg .doneCrit0,#mermaid-svg-4vuN1MjXmBoSUadg .doneCrit1,#mermaid-svg-4vuN1MjXmBoSUadg .doneCrit2,#mermaid-svg-4vuN1MjXmBoSUadg .doneCrit3{stroke:#f88;fill:#d3d3d3;stroke-width:2;cursor:pointer;shape-rendering:crispEdges}#mermaid-svg-4vuN1MjXmBoSUadg .milestone{transform:rotate(45deg) scale(0.8, 0.8)}#mermaid-svg-4vuN1MjXmBoSUadg .milestoneText{font-style:italic}#mermaid-svg-4vuN1MjXmBoSUadg .doneCritText0,#mermaid-svg-4vuN1MjXmBoSUadg .doneCritText1,#mermaid-svg-4vuN1MjXmBoSUadg .doneCritText2,#mermaid-svg-4vuN1MjXmBoSUadg .doneCritText3{fill:#000 !important}#mermaid-svg-4vuN1MjXmBoSUadg .activeCritText0,#mermaid-svg-4vuN1MjXmBoSUadg .activeCritText1,#mermaid-svg-4vuN1MjXmBoSUadg .activeCritText2,#mermaid-svg-4vuN1MjXmBoSUadg .activeCritText3{fill:#000 !important}#mermaid-svg-4vuN1MjXmBoSUadg .titleText{text-anchor:middle;font-size:18px;fill:#000;font-family:'trebuchet ms', verdana, arial;font-family:var(--mermaid-font-family)}#mermaid-svg-4vuN1MjXmBoSUadg g.classGroup text{fill:#9370db;stroke:none;font-family:'trebuchet ms', verdana, arial;font-family:var(--mermaid-font-family);font-size:10px}#mermaid-svg-4vuN1MjXmBoSUadg g.classGroup text .title{font-weight:bolder}#mermaid-svg-4vuN1MjXmBoSUadg g.clickable{cursor:pointer}#mermaid-svg-4vuN1MjXmBoSUadg g.classGroup rect{fill:#ECECFF;stroke:#9370db}#mermaid-svg-4vuN1MjXmBoSUadg g.classGroup line{stroke:#9370db;stroke-width:1}#mermaid-svg-4vuN1MjXmBoSUadg .classLabel .box{stroke:none;stroke-width:0;fill:#ECECFF;opacity:0.5}#mermaid-svg-4vuN1MjXmBoSUadg .classLabel .label{fill:#9370db;font-size:10px}#mermaid-svg-4vuN1MjXmBoSUadg .relation{stroke:#9370db;stroke-width:1;fill:none}#mermaid-svg-4vuN1MjXmBoSUadg .dashed-line{stroke-dasharray:3}#mermaid-svg-4vuN1MjXmBoSUadg #compositionStart{fill:#9370db;stroke:#9370db;stroke-width:1}#mermaid-svg-4vuN1MjXmBoSUadg #compositionEnd{fill:#9370db;stroke:#9370db;stroke-width:1}#mermaid-svg-4vuN1MjXmBoSUadg #aggregationStart{fill:#ECECFF;stroke:#9370db;stroke-width:1}#mermaid-svg-4vuN1MjXmBoSUadg #aggregationEnd{fill:#ECECFF;stroke:#9370db;stroke-width:1}#mermaid-svg-4vuN1MjXmBoSUadg #dependencyStart{fill:#9370db;stroke:#9370db;stroke-width:1}#mermaid-svg-4vuN1MjXmBoSUadg #dependencyEnd{fill:#9370db;stroke:#9370db;stroke-width:1}#mermaid-svg-4vuN1MjXmBoSUadg #extensionStart{fill:#9370db;stroke:#9370db;stroke-width:1}#mermaid-svg-4vuN1MjXmBoSUadg #extensionEnd{fill:#9370db;stroke:#9370db;stroke-width:1}#mermaid-svg-4vuN1MjXmBoSUadg .commit-id,#mermaid-svg-4vuN1MjXmBoSUadg .commit-msg,#mermaid-svg-4vuN1MjXmBoSUadg .branch-label{fill:lightgrey;color:lightgrey;font-family:'trebuchet ms', verdana, arial;font-family:var(--mermaid-font-family)}#mermaid-svg-4vuN1MjXmBoSUadg .pieTitleText{text-anchor:middle;font-size:25px;fill:#000;font-family:'trebuchet ms', verdana, arial;font-family:var(--mermaid-font-family)}#mermaid-svg-4vuN1MjXmBoSUadg .slice{font-family:'trebuchet ms', verdana, arial;font-family:var(--mermaid-font-family)}#mermaid-svg-4vuN1MjXmBoSUadg g.stateGroup text{fill:#9370db;stroke:none;font-size:10px;font-family:'trebuchet ms', verdana, arial;font-family:var(--mermaid-font-family)}#mermaid-svg-4vuN1MjXmBoSUadg g.stateGroup text{fill:#9370db;fill:#333;stroke:none;font-size:10px}#mermaid-svg-4vuN1MjXmBoSUadg g.statediagram-cluster .cluster-label text{fill:#333}#mermaid-svg-4vuN1MjXmBoSUadg g.stateGroup .state-title{font-weight:bolder;fill:#000}#mermaid-svg-4vuN1MjXmBoSUadg g.stateGroup rect{fill:#ECECFF;stroke:#9370db}#mermaid-svg-4vuN1MjXmBoSUadg g.stateGroup line{stroke:#9370db;stroke-width:1}#mermaid-svg-4vuN1MjXmBoSUadg .transition{stroke:#9370db;stroke-width:1;fill:none}#mermaid-svg-4vuN1MjXmBoSUadg .stateGroup .composit{fill:white;border-bottom:1px}#mermaid-svg-4vuN1MjXmBoSUadg .stateGroup .alt-composit{fill:#e0e0e0;border-bottom:1px}#mermaid-svg-4vuN1MjXmBoSUadg .state-note{stroke:#aa3;fill:#fff5ad}#mermaid-svg-4vuN1MjXmBoSUadg .state-note text{fill:black;stroke:none;font-size:10px}#mermaid-svg-4vuN1MjXmBoSUadg .stateLabel .box{stroke:none;stroke-width:0;fill:#ECECFF;opacity:0.7}#mermaid-svg-4vuN1MjXmBoSUadg .edgeLabel text{fill:#333}#mermaid-svg-4vuN1MjXmBoSUadg .stateLabel text{fill:#000;font-size:10px;font-weight:bold;font-family:'trebuchet ms', verdana, arial;font-family:var(--mermaid-font-family)}#mermaid-svg-4vuN1MjXmBoSUadg .node circle.state-start{fill:black;stroke:black}#mermaid-svg-4vuN1MjXmBoSUadg .node circle.state-end{fill:black;stroke:white;stroke-width:1.5}#mermaid-svg-4vuN1MjXmBoSUadg #statediagram-barbEnd{fill:#9370db}#mermaid-svg-4vuN1MjXmBoSUadg .statediagram-cluster rect{fill:#ECECFF;stroke:#9370db;stroke-width:1px}#mermaid-svg-4vuN1MjXmBoSUadg .statediagram-cluster rect.outer{rx:5px;ry:5px}#mermaid-svg-4vuN1MjXmBoSUadg .statediagram-state .divider{stroke:#9370db}#mermaid-svg-4vuN1MjXmBoSUadg .statediagram-state .title-state{rx:5px;ry:5px}#mermaid-svg-4vuN1MjXmBoSUadg .statediagram-cluster.statediagram-cluster .inner{fill:white}#mermaid-svg-4vuN1MjXmBoSUadg .statediagram-cluster.statediagram-cluster-alt .inner{fill:#e0e0e0}#mermaid-svg-4vuN1MjXmBoSUadg .statediagram-cluster .inner{rx:0;ry:0}#mermaid-svg-4vuN1MjXmBoSUadg .statediagram-state rect.basic{rx:5px;ry:5px}#mermaid-svg-4vuN1MjXmBoSUadg .statediagram-state rect.divider{stroke-dasharray:10,10;fill:#efefef}#mermaid-svg-4vuN1MjXmBoSUadg .note-edge{stroke-dasharray:5}#mermaid-svg-4vuN1MjXmBoSUadg .statediagram-note rect{fill:#fff5ad;stroke:#aa3;stroke-width:1px;rx:0;ry:0}:root{--mermaid-font-family: '"trebuchet ms", verdana, arial';--mermaid-font-family: "Comic Sans MS", "Comic Sans", cursive}#mermaid-svg-4vuN1MjXmBoSUadg .error-icon{fill:#522}#mermaid-svg-4vuN1MjXmBoSUadg .error-text{fill:#522;stroke:#522}#mermaid-svg-4vuN1MjXmBoSUadg .edge-thickness-normal{stroke-width:2px}#mermaid-svg-4vuN1MjXmBoSUadg .edge-thickness-thick{stroke-width:3.5px}#mermaid-svg-4vuN1MjXmBoSUadg .edge-pattern-solid{stroke-dasharray:0}#mermaid-svg-4vuN1MjXmBoSUadg .edge-pattern-dashed{stroke-dasharray:3}#mermaid-svg-4vuN1MjXmBoSUadg .edge-pattern-dotted{stroke-dasharray:2}#mermaid-svg-4vuN1MjXmBoSUadg .marker{fill:#333}#mermaid-svg-4vuN1MjXmBoSUadg .marker.cross{stroke:#333}:root { --mermaid-font-family: "trebuchet ms", verdana, arial;}#mermaid-svg-4vuN1MjXmBoSUadg {color: rgba(0, 0, 0, 0.75);font: ;}

终止条件

其他情况

Soft Thresholding

python 实现：
in ista.py

from utils import sh
import numpy as npclass ISTA:_w_d = np.asmatrix([])_lambda = 0_l = 0_alpha = 0'''init(wd, lam, l)wd      字典      wd      lam     系数      lambda  l       常数      l       '''def __init__(self, wd, lam, lip, alpha):self._w_d = np.asmatrix(wd)self._lambda = lameig, eig_vs = np.linalg.eig(np.dot(wd.T, wd))r = np.max(eig)if r > lip:print('L is not suitable where l = %f but eig_max(wd) = %f' % (lip, r))print('now L is changed to eig_max(wd) + 1 which larger than eig_max(wd)')lip = r + 1del eig, eig_vsself._l = lipself._alpha = alphaself.we = self._w_d.T / self._lself.beta = self._alpha / self._l'''y,e = do_train(x, eps, itm)x       原向量         x                eps     误差限         epsilon          itm     最大迭代次数    iterate_times_max'''def _work(self, x, eps):z = np.zeros((self._w_d.shape[1], 1))it = 0err_t = []while True:it += 1t = np.dot(self._w_d, z) - xz1 = sh(z - self.we * t, self.beta)err = np.linalg.norm(z - z1)err_t.append(err)if err < eps:breakz = z1return z, it, np.asarray(err_t)def do_work(self, args):assert args.__len__() == 2return self._work(args[0], args[1])

in utils.py

import numpy as npdef sh(x, theta):return np.multiply(np.sign(x), np.maximum(np.abs(x) - theta, 0))def model_work(model, *args):return model.do_work(args);

in entrance.py

from utils import model_train, model_work
from ista import ISTA
import numpy as np
import tensorflow as tf
from scipy.linalg import orthdef get_wd(m, n):psi = np.eye(m)phi = np.random.randn(n, m)phi = np.transpose(orth(np.transpose(phi)))return np.dot(phi, psi)def get_x_z(m, k, wd):z = np.zeros([m, 1])ind = np.random.choice(a=m, size=k, replace=False, p=None)z[ind, 0] = np.random.randn(1, k) * 10return np.dot(wd, z), zdef call_ista(wd, lam, lip, alpha, x, eps):model = ISTA(wd, lam, lip, alpha)return model_work(model, x, eps)def test_ista(m, n, k, wd, x, z, lam, lip, alpha, eps):zr, it, err = call_ista(wd, lam, lip, alpha, x, eps)with open('ista.txt', 'w') as f:for i in range(m):f.write('%03.3f %03.3f\n' % (z[i, 0], zr[i, 0]))print(it)print(err)def train_lista(m, n, k, wd, lam, L, alpha, eps, times):# model = Lista(...)for i in range(times):x, z = get_x_z(m, k, wd)# train_model(model,x,z,T)passdef main():m, n, k = 1024, 256, 8wd = get_wd(m, n)x, z = get_x_z(m, k, wd)L = 2lam = 0.1alpha = 0.1eps = 0.001test_ista(m, n, k, wd, x, z, lam, L, alpha, eps)if __name__ == '__main__':main()

matlab代码

% CSDN不支持matlab代码，这里使用lua近似替代，可能有显示问题
% author sunwei
% at 2020.10.11
function z = sh(beta, theta)
% sh为软阈值函数简写
% 传入一个向量beta，一个阈值theta
% z 为软阈值变换的结果
z = sign(beta).*max(abs(beta)-theta,0);
endfunction [x,z] = get_x_z(m,k,wd)
% 随机信号生成
% z为`稀疏信号`
% x为`原始信号`
% 为方便验证，实际上原始信号通过稀疏信号生成
z = zeros(m,1);
ind = randperm(m);
ind = ind(1:k);
z(ind) = randn(1,k) * 10;
x = wd*z;
endfunction [zr, it] = ista(wd, L, alpha, x, eps)
% wd 字典
% alpha sh参数的一部分
% L 李普希斯常数(一个大于wd'*wd)的常数
% x 得到的信号
% eps 误差[~,n] = size(wd);z = zeros(n,1);it = 0;while truet = z - 1/L*wd'*(wd*z-x);zr = sh(t, alpha/L);it = it + 1;if((z-zr)'*(z-zr) < eps) break; endz = zr;end
endfunction test_ista()m = 1024; n = 256; k = 6;wd = randn(n,m);wd = orth(wd')';[x, z] = get_x_z(m,k,wd);eig_val = eig(wd'*wd);L = max(eig_val) + 1;lam = 0.1;alpha = 0.1;eps = 0.001;[zr, it] = ista(wd, L, alpha, x, eps);t = 1:m;plot(t,z,t,zr)fprintf('it = %d, err = %d', it, (z-zr)'*(z-zr));
end

补充：

下一篇：LISTA 基于学习的迭代压缩软阈值算法(Paper: Learning Fast Approximations of Sparse Coding)
关于稀疏表示可以看一看Coursera上的https://www.coursera.org/learn/image-processing，Duke大学的Professor还是厉害的。
python fprintf('it = %d, err = %d', it, (z-zr)'*(z-zr));一句有问题，err应该是稀疏性提供的误差+保真度提供的误差