A using case: iterative dominant eigensolver

[GiggleLiu]

using Zygote, NiSparseArrays
using SparseArrays, LinearAlgebra

function iterative_ground_state_solver(A, x, target; niter=100)
    for i=1:niter
        x = A * x
        x /= norm(x)
    end
    return abs(x' * target)
end

A = sprand(100, 100, 0.1)
x = randn(100)
target = randn(100)
ga, gx, gt = Zygote.gradient(iterative_ground_state_solver, A, x, target)

Here, I first compute the dominant eigenvalue of A, and use the fidelity between this dominant eigenvalue and target vector as the loss.

[jieli-matrix]

这是个很好的例子，但我感觉可能没有完全解决我的问题：
我不太懂现在这个高阶算子的微分目标是什么r？
以幂法为例，输入是随机的向量x和给定矩阵A，输出是矩阵A的模最大特征向量；
里面涉及到的操作是 矩阵-向量乘 , norm(x)计算 和 点除；是要对每一步都微分？
但目前rrule只导出了矩阵-向量乘的规则
我要写的function iterative_ground_state_solver 是NiLang函数，还是普通函数？
之前的rrule怎么用上？
后续的svd_lowrank会涉及到QR与svd操作，这些都是没有实现的，它们的rrule怎么处理？
cc: @johnnychen94

[GiggleLiu]

这个高阶算子的微分目标是什么r？

没听懂，你说目标函数？是fidelity，即和target的相似度。

里面涉及到的操作是 矩阵-向量乘 , norm(x)计算 和 点除；是要对每一步都微分？
后续的svd_lowrank会涉及到QR与svd操作，这些都是没有实现的，它们的rrule怎么处理？

svd, norm, 点除这些的AD rule 都是ChainRules里面现成的。每一步都需要微分，但是Zygote自动帮你串联起来了。

我要写的function iterative_ground_state_solver 是NiLang函数，还是普通函数？

普通函数

之前的rrule怎么用上？

你是说你写的那些rrule吗？using 了就可以了。

[jieli-matrix]

此外，幂法只能求解模最大特征值/特征向量；如果想要求解学长例子里的最小特征值对应的特征向量，需要采用反幂法，对inv(A)做幂法；这就会涉及到稀疏矩阵的LU分解。ChainRules会包含lu和qr的AD rule吗？还是说这两个操作不需要微分？
以及如果求解的A如学长生成的几乎近似三对角矩阵，为什么不用对称QR方法去求解特征值呢？

[GiggleLiu]

此外，幂法只能求解模最大特征值/特征向量；如果想要求解学长例子里的最小特征值对应的特征向量，需要采用反幂法，对inv(A)做幂法；

可以求最小的，只要做个shift就行，类似加个identity矩阵，在加个负号，很容易把最小eigenvalue变成最大的，inv不是必要的。

ChainRules会包含lu和qr的AD rule吗？

具体我记不清，但是即使没有的话，另一个仓库抄过来就行：
https://github.com/GiggleLiu/BackwardsLinalg.jl
这些方法都是dense矩阵基础上推导的～不一定可以很容易拓展到稀疏矩阵。

[jieli-matrix]

那样就没问题，最小最大的幂法都可以写，我去试一下

但我刚刚查了...chainrules没有qr的AD rule...
https://github.com/JuliaDiff/ChainRules.jl/blob/master/src/rulesets/LinearAlgebra/factorization.jl
学长给BackwardsLinalg.jl加下稀疏矩阵qr的支持，然后我把svd_lowrank扩展到稀疏矩阵应该没问题

[GiggleLiu]

学长给BackwardsLinalg.jl加下稀疏矩阵qr的支持

要不你试着自己给ChainRules提个PR? 正好是个锻炼的机会～

[jieli-matrix]

如果有这样的能力，我自然是愿意去做的；但我觉得恐怕时间来不及，这方面微分规则推导我也不熟悉；而且项目规划里是在low-level的稀疏矩阵乘法上实现svd_lowrank，所以我当时以为qr和svd应该是有配套的设施不需要我再去做build-in的低效rules实现？

[GiggleLiu]

行吧，明明已经把backward rule实现都贴给你了，我来提PR吧～

[johnnychen94]

稀疏和非稀疏的情况在实现上应该差别还是蛮大的. 再者试图构造稀疏矩阵QR的求导真的有意义么？QR分解的结果也与稀疏矩阵的0值有关，所以如果要正确求导的话，应该是要考虑0值的导数。这种情况下稀疏的qr分解求导其实一定程度上已经退化成了稠密的qr分解求导。

[GiggleLiu]

@johnnychen94 这里我们其实讨论了两个问题，实现low rank svd应该只涉及dense矩阵的qr：https://github.com/jieli-matrix/svd_lowrank/blob/610cf33ac59ad9235fa5cc57f166a55ed01d9e1e/pytorch_implement/get_approximate_basis.py#L50

学妹之前的顾虑是求稀疏矩阵逆过程中需要稀疏矩阵QR，稀疏矩阵的逆是个很不好做的问题，因为结果无法直接表示为稀疏矩阵，或许可以在krylov空间做，或者一些对接近对角的稀疏矩阵作近似。但这些不重要，因为可以绕开，很少有人会直接求稀疏矩阵逆。

[jieli-matrix]

我今早想了一下，low rank svd确实只涉及到dense矩阵的qr，所以目前是需要封装成rrule并通过CRTU的测试
稀疏矩阵求逆是用在反幂法上，也不是直接求逆，做的是sparse lu，学长说的那种方法求最小特征值是可以，但收敛速度可能会慢，因为shift需要多次尝试

[GiggleLiu]

我和久宁聊了下，感觉没有demo别人用这个库可能有点困难，
要不你把我写的这个例子当作using case加进README里面当作小的demo吧，也许没有必要写一个很完整的demo。

[jieli-matrix]

Closed since I open it in #33