Source code for cil.optimisation.operators.DiagonalOperator
# Copyright 2020 United Kingdom Research and Innovation
# Copyright 2020 The University of Manchester
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
# http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.
#
# Authors:
# CIL Developers, listed at: https://github.com/TomographicImaging/CIL/blob/master/NOTICE.txt
import numpy as np
from cil.framework import ImageData
from cil.optimisation.operators import LinearOperator
[docs]
class DiagonalOperator(LinearOperator):
r"""DiagonalOperator
Performs an element-wise multiplication, i.e., `Hadamard Product <https://en.wikipedia.org/wiki/Hadamard_product_(matrices)#:~:text=In%20mathematics%2C%20the%20Hadamard%20product,elements%20i%2C%20j%20of%20the>`_
of a :class:`DataContainer` `x` and :class:`DataContainer` `diagonal`, `d` .
.. math:: (D\circ x) = \sum_{i,j}^{M,N} D_{i,j} x_{i, j}
In matrix-vector interpretation, if `D` is a :math:`M\times N` dense matrix and is flattened, we have a :math:`M*N \times M*N` vector.
A sparse diagonal matrix, i.e., :class:`DigaonalOperator` can be created if we add the vector above to the main diagonal.
If the :class:`DataContainer` `x` is also flattened, we have a :math:`M*N` vector.
Now, matrix-vector multiplcation is allowed and results to a :math:`(M*N,1)` vector. After reshaping we recover a :math:`M\times N` :class:`DataContainer`.
Parameters
----------
diagonal : DataContainer
DataContainer with the same dimensions as the data to be operated on.
domain_geometry : ImageGeometry
Specifies the geometry of the operator domain. If 'None' will use the diagonal geometry directly. default=None .
"""
def __init__(self, diagonal, domain_geometry=None):
if domain_geometry is None:
domain_geometry = diagonal.geometry.copy()
super(DiagonalOperator, self).__init__(domain_geometry=domain_geometry,
range_geometry=domain_geometry)
self.diagonal = diagonal
[docs]
def direct(self,x,out=None):
"Returns :math:`D\circ x` "
if out is None:
return self.diagonal * x
else:
self.diagonal.multiply(x,out=out)
return out
[docs]
def adjoint(self,x, out=None):
"Returns :math:`D^*\circ x` "
return self.diagonal.conjugate().multiply(x,out=out)
[docs]
def calculate_norm(self, **kwargs):
r""" Returns the operator norm of DiagonalOperator which is the :math:`\infty` norm of `diagonal`
.. math:: \|D\|_{\infty} = \max_{i}\{|D_{i}|\}
"""
return self.diagonal.abs().max()