some problem about the discrete of the first derivative operator

4 次查看（过去 30 天）
bingbing
bingbing 2024-3-5
评论： bingbing 2024-3-5
I am reading a paper(Parameter Choice Strategies for Multipenalty Regularization Massimo Fornasier, Valeriya Naumova, and Sergei V. Pereverzyev SIAM Journal on Numerical Analysis 2014 52:4, 1770-1794)
In this paper first derivative operator is a matrix like the pics below
is this right?
because most of the first-order derivative discrete operator matrices I've seen are like this
from(Discrete inverse problems: insight and algorithms,Hansen)
Could someone please explain this to me?
And how to calculate the matrix B in the first pics with matlab?The B matrix I obtained using the "sqrtm" command is a complex matrix, is this right?

请先登录，再回答此问题。

采纳的回答

John D'Errico
John D'Errico 2024-3-5
编辑：John D'Errico 2024-3-5
Ran in:
It is difficult to know, not without reading the article.
Is the D matrix as shown just a typo that got by the proofreaders? Very possibly, but as I show below, it may be irrelevant. Assuming the signs on alternate rows are flipped, this would generate alternating signs on the derivative approximations, which seems a bit strange. I don't see any valid reason for it, but then I don't have the paper in front of me either.
n = 10;
L = -diag(ones(n,1)) + diag(ones(n-1,1),1);
L(end,:) = []
L = 9×10
-1 1 0 0 0 0 0 0 0 0 0 -1 1 0 0 0 0 0 0 0 0 0 -1 1 0 0 0 0 0 0 0 0 0 -1 1 0 0 0 0 0 0 0 0 0 -1 1 0 0 0 0 0 0 0 0 0 -1 1 0 0 0 0 0 0 0 0 0 -1 1 0 0 0 0 0 0 0 0 0 -1 1 0 0 0 0 0 0 0 0 0 -1 1
But we can create the matrix D which alternates as their matrix seems to do simply from L, as you wrote it:
T = diag(mod(0:n-2,2)*2-1)
T = 9×9
-1 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 -1 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 -1 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 -1 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 -1
D = T*L
D = 9×10
1 -1 0 0 0 0 0 0 0 0 0 -1 1 0 0 0 0 0 0 0 0 0 1 -1 0 0 0 0 0 0 0 0 0 -1 1 0 0 0 0 0 0 0 0 0 1 -1 0 0 0 0 0 0 0 0 0 -1 1 0 0 0 0 0 0 0 0 0 1 -1 0 0 0 0 0 0 0 0 0 -1 1 0 0 0 0 0 0 0 0 0 1 -1
In the end, does it matter? No, at least not here. This is because D'*D == L'*L. Yes. The prodiuct will be identically the same.
Look carefully at how I defined D = T*L. Expand the product D'*D, and you will see that
D'*D = (T*L)' * (T*L) = L'*T'*T*L = L'*L
since T is a diagonal matrix, composed on 1 and -1 on the diagonal. When you form T'*T, that turns into the identity matrix.
As such, while I have no clue if they intended to write that matrix for D properly or it was just a typo, if the only thing they will do with it is to create a penalizing operator, then it is irrelevant.
Now, is it correct that the B matrix is complex? Well, yes, and no. I want to say that is a complex question, but then I expect at least some of my readers will start to groan. The matrix
DprimeD = D'*D
DprimeD = 10×10
1 -1 0 0 0 0 0 0 0 0 -1 2 -1 0 0 0 0 0 0 0 0 -1 2 -1 0 0 0 0 0 0 0 0 -1 2 -1 0 0 0 0 0 0 0 0 -1 2 -1 0 0 0 0 0 0 0 0 -1 2 -1 0 0 0 0 0 0 0 0 -1 2 -1 0 0 0 0 0 0 0 0 -1 2 -1 0 0 0 0 0 0 0 0 -1 2 -1 0 0 0 0 0 0 0 0 -1 1
is rank deficient. It MUST be so. It will have one zero eigenvalue, since the matrix D'*D
D'*D
ans = 10×10
1 -1 0 0 0 0 0 0 0 0 -1 2 -1 0 0 0 0 0 0 0 0 -1 2 -1 0 0 0 0 0 0 0 0 -1 2 -1 0 0 0 0 0 0 0 0 -1 2 -1 0 0 0 0 0 0 0 0 -1 2 -1 0 0 0 0 0 0 0 0 -1 2 -1 0 0 0 0 0 0 0 0 -1 2 -1 0 0 0 0 0 0 0 0 -1 2 -1 0 0 0 0 0 0 0 0 -1 1
will be an nxn matrix, but it was created as a product of matrices nx(n-1) matrices, each of which will have rank exactly n-1. So we see that
rank(DprimeD)
ans = 9
will also be n-1, as it was in my example.
Now the problem arises that sqrtm gets a little lost. Effectively, sqrtm does not know this matrix is actually rank deficient. It just does the best it can in double precision arithmetic, and along the way, it generates some TINY imaginary parts. The result is an array with infinitessimal imaginary parts.
sqrtm(DprimeD)
ans =
0.8436 + 0.0000i -0.5146 + 0.0000i -0.1267 + 0.0000i -0.0622 + 0.0000i -0.0389 + 0.0000i -0.0279 + 0.0000i -0.0220 + 0.0000i -0.0186 + 0.0000i -0.0168 + 0.0000i -0.0159 + 0.0000i -0.5146 + 0.0000i 1.2315 + 0.0000i -0.4501 + 0.0000i -0.1034 + 0.0000i -0.0512 + 0.0000i -0.0330 + 0.0000i -0.0245 + 0.0000i -0.0201 + 0.0000i -0.0178 + 0.0000i -0.0168 + 0.0000i -0.1267 + 0.0000i -0.4501 + 0.0000i 1.2548 + 0.0000i -0.4391 + 0.0000i -0.0975 + 0.0000i -0.0478 + 0.0000i -0.0311 + 0.0000i -0.0237 + 0.0000i -0.0201 + 0.0000i -0.0186 + 0.0000i -0.0622 + 0.0000i -0.1034 + 0.0000i -0.4391 + 0.0000i 1.2607 + 0.0000i -0.4358 + 0.0000i -0.0956 + 0.0000i -0.0470 + 0.0000i -0.0311 + 0.0000i -0.0245 + 0.0000i -0.0220 + 0.0000i -0.0389 + 0.0000i -0.0512 + 0.0000i -0.0975 + 0.0000i -0.4358 + 0.0000i 1.2626 + 0.0000i -0.4349 + 0.0000i -0.0956 + 0.0000i -0.0478 + 0.0000i -0.0330 + 0.0000i -0.0279 + 0.0000i -0.0279 + 0.0000i -0.0330 + 0.0000i -0.0478 + 0.0000i -0.0956 + 0.0000i -0.4349 + 0.0000i 1.2626 + 0.0000i -0.4358 + 0.0000i -0.0975 + 0.0000i -0.0512 + 0.0000i -0.0389 + 0.0000i -0.0220 + 0.0000i -0.0245 + 0.0000i -0.0311 + 0.0000i -0.0470 + 0.0000i -0.0956 + 0.0000i -0.4358 + 0.0000i 1.2607 + 0.0000i -0.4391 + 0.0000i -0.1034 + 0.0000i -0.0622 + 0.0000i -0.0186 + 0.0000i -0.0201 + 0.0000i -0.0237 + 0.0000i -0.0311 + 0.0000i -0.0478 + 0.0000i -0.0975 + 0.0000i -0.4391 + 0.0000i 1.2548 + 0.0000i -0.4501 + 0.0000i -0.1267 + 0.0000i -0.0168 + 0.0000i -0.0178 + 0.0000i -0.0201 + 0.0000i -0.0245 + 0.0000i -0.0330 + 0.0000i -0.0512 + 0.0000i -0.1034 + 0.0000i -0.4501 + 0.0000i 1.2315 + 0.0000i -0.5146 + 0.0000i -0.0159 + 0.0000i -0.0168 + 0.0000i -0.0186 + 0.0000i -0.0220 + 0.0000i -0.0279 + 0.0000i -0.0389 + 0.0000i -0.0622 + 0.0000i -0.1267 + 0.0000i -0.5146 + 0.0000i 0.8436 + 0.0000i
Really, this is just due to tiny rounding errors in the least significant bits. You can choose to discard the imaginary parts, since we know they are just floating point trash. Well, I know that. I can't guess what you know, but then, I'm the one who is writing this response. ;-)
B = real(sqrtm(DprimeD))
B = 10×10
0.8436 -0.5146 -0.1267 -0.0622 -0.0389 -0.0279 -0.0220 -0.0186 -0.0168 -0.0159 -0.5146 1.2315 -0.4501 -0.1034 -0.0512 -0.0330 -0.0245 -0.0201 -0.0178 -0.0168 -0.1267 -0.4501 1.2548 -0.4391 -0.0975 -0.0478 -0.0311 -0.0237 -0.0201 -0.0186 -0.0622 -0.1034 -0.4391 1.2607 -0.4358 -0.0956 -0.0470 -0.0311 -0.0245 -0.0220 -0.0389 -0.0512 -0.0975 -0.4358 1.2626 -0.4349 -0.0956 -0.0478 -0.0330 -0.0279 -0.0279 -0.0330 -0.0478 -0.0956 -0.4349 1.2626 -0.4358 -0.0975 -0.0512 -0.0389 -0.0220 -0.0245 -0.0311 -0.0470 -0.0956 -0.4358 1.2607 -0.4391 -0.1034 -0.0622 -0.0186 -0.0201 -0.0237 -0.0311 -0.0478 -0.0975 -0.4391 1.2548 -0.4501 -0.1267 -0.0168 -0.0178 -0.0201 -0.0245 -0.0330 -0.0512 -0.1034 -0.4501 1.2315 -0.5146 -0.0159 -0.0168 -0.0186 -0.0220 -0.0279 -0.0389 -0.0622 -0.1267 -0.5146 0.8436
norm(imag(sqrtm(DprimeD)))
ans = 2.2726e-08
And, working in double precision arithmetic, I would expect that norm to be on the order of sqrt(eps). As you can see, it was very much so.
sqrt(eps)
ans = 1.4901e-08

更多回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 Multidimensional Arrays 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by