diff --git a/pyrobolearn/backends/torch_backend.py b/pyrobolearn/backends/torch_backend.py index 8d43fa7..7ee1b0b 100644 --- a/pyrobolearn/backends/torch_backend.py +++ b/pyrobolearn/backends/torch_backend.py @@ -97,6 +97,8 @@ def concatenate(data, axis=0, out=None): # missing torch.dstack # missing torch.vstack # missing torch.hstack +# np.ndim --> torch.dim() +# np.reshape --> torch.view() # Tests diff --git a/pyrobolearn/losses/__init__.py b/pyrobolearn/losses/__init__.py index c21dfc1..5b893a6 100644 --- a/pyrobolearn/losses/__init__.py +++ b/pyrobolearn/losses/__init__.py @@ -1,3 +1,4 @@ # import losses from .loss import * +from .losses import * diff --git a/pyrobolearn/losses/loss.py b/pyrobolearn/losses/loss.py index 2bd3d1f..4250ad7 100644 --- a/pyrobolearn/losses/loss.py +++ b/pyrobolearn/losses/loss.py @@ -1,13 +1,13 @@ #!/usr/bin/env python -"""Defines the common loss functions that are used by the learning algorithm / optimizer. +"""Defines the abstract loss class that is used by the learning algorithm / optimizer. Losses are evaluated on model parameters, data batches, and / or storages. -TODO: loss vs cost vs reward -costs and rewards are defined possibly for each time steps. +Note that rewards / costs and losses are different concepts. Losses are minimized with respect to parameters to +optimize them, while rewards / costs depends on the state(s) and action(s) and are returned by the environment. """ -from abc import ABCMeta, abstractmethod +from abc import ABCMeta import operator import copy import collections @@ -283,229 +283,6 @@ class Loss(object): return loss -class FixedLoss(Loss): - r"""Fixed Loss - - """ - def __init__(self, value): - super(FixedLoss, self).__init__() - self.value = value - - def compute(self, batch): - return self.value * batch - - -class L2Loss(Loss): - r"""L2 Loss - - Compute the L2 loss given by: :math:`1/2 * (y_{target} - y_{predict})^2` - """ - def __init__(self, target, approximator): - super(L2Loss, self).__init__() - self.target = target - self.approximator = approximator - - def compute(self, batch): - # based on approximator check what we need - return 0.5 * (self.target(arg) - self.approximator(arg)).pow(2).mean() - - -class ValueLoss(Loss): - r"""L2 loss for values - """ - def __init__(self): - super(ValueLoss, self).__init__() - - def compute(self, batch): - returns = batch['returns'] - values = batch.current['values'] - return 0.5 * (returns - values).pow(2).mean() - - -class PGLoss(Loss): - r"""Policy Gradient Loss - - Compute the policy gradient loss which is maximized and given by: - - .. math:: L^{PG} = \mathbb{E}[ \log \pi_{\theta}(a_t | s_t) \psi_t ] - - where :math:`\psi_t` is the associated return estimator, which can be for instance, the total reward estimator - :math:`\psi_t = R(\tau)` (where :math:`\tau` represents the whole trajectory), the state action value estimator - :math:`\psi_t = Q(s_t, a_t)`, or the advantage estimator :math:`\psi_t = A_t = Q(s_t, a_t) - V(s_t)`. Other - estimators are also possible. - - The gradient with respect to the parameters :math:`\theta` is then given by: - - .. math:: g = \mathbb{E}[ \nabla_\theta \log \pi_{\theta}(a_t | s_t) - - References: - [1] "Proximal Policy Optimization Algorithms", Schulman et al., 2017 - [2] "High-Dimensional Continuous Control using Generalized Advantage Estimation", Schulman et al., 2016 - """ - - def __init__(self): - super(PGLoss, self).__init__() - - def compute(self, batch): - log_curr_pi = batch.current['action_distributions'] - log_curr_pi = log_curr_pi.log_probs(batch.current['actions']) - estimator = batch['estimator'] - loss = torch.exp(log_curr_pi) * estimator - return -loss.mean() - - def latex(self): - return "\\mathbb{E}[ r_t(\\theta) A_t ]" - - -class CPILoss(Loss): - r"""CPI Loss - - Conservative Policy Iteration objective which is maximized and defined in [1]: - - .. math:: L^{CPI}(\theta) = \mathbb{E}[ r_t(\theta) A_t ] - - where the expectation is taken over a finite batch of samples, :math:`A_t` is an estimator of the advantage fct at - time step :math:`t`, :math:`r_t(\theta)` is the probability ratio given by - :math:`r_t(\theta) = \frac{ \pi_{\theta}(a_t|s_t) }{ \pi_{\theta_{old}}(a_t|s_t) }`. - - References: - [1] "Approximately optimal approximate reinforcement learning", Kakade et al., 2002 - [2] "Proximal Policy Optimization Algorithms", Schulman et al., 2017 - """ - - def __init__(self): - """ - Initialize the CPI Loss. - """ - super(CPILoss, self).__init__() - - def compute(self, batch): # policy_distribution, old_policy_distribution, estimator): - # ratio = policy_distribution / old_policy_distribution - log_curr_pi = batch.current['action_distributions'] - log_curr_pi = log_curr_pi.log_probs(batch.current['actions']) - log_prev_pi = batch['action_distributions'] - log_prev_pi = log_prev_pi.log_probs(batch['actions']) - - ratio = torch.exp(log_curr_pi - log_prev_pi) - estimator = batch['estimator'] - - loss = ratio * estimator - return -loss.mean() - - def latex(self): - return "\\mathbb{E}[ r_t(\\theta) A_t ]" - - -class CLIPLoss(Loss): - r"""CLIP Loss - - Loss defined in [1] which is maximized and given by: - - .. math:: L^{CLIP}(\theta) = \mathbb{E}[ \min(r_t(\theta) A_t, clip(r_t(\theta), 1-\epsilon, 1+\epsilon) A_t) ] - - where the expectation is taken over a finite batch of samples, :math:`A_t` is an estimator of the advantage fct at - time step :math:`t`, :math:`r_t(\theta)` is the probability ratio given by - :math:`r_t(\theta) = \frac{ \pi_{\theta}(a_t|s_t) }{ \pi_{\theta_{old}}(a_t|s_t) }`. - - References: - [1] "Proximal Policy Optimization Algorithms", Schulman et al., 2017 - """ - - def __init__(self, clip=0.2): - """ - Initialize the loss. - - Args: - epsilon (float): clip parameter - """ - super(CLIPLoss, self).__init__() - self.eps = clip - - def compute(self, batch): # , policy_distribution, old_policy_distribution, estimator): - log_curr_pi = batch.current['action_distributions'] - log_curr_pi = log_curr_pi.log_probs(batch.current['actions']) - log_prev_pi = batch['action_distributions'] - log_prev_pi = log_prev_pi.log_probs(batch['actions']) - - ratio = torch.exp(log_curr_pi - log_prev_pi) - estimator = batch['estimator'] - - loss = torch.min(ratio * estimator, torch.clamp(ratio, 1.0-self.eps, 1.0+self.eps) * estimator) - return -loss.mean() - - def latex(self): - return "\\mathbb{E}[ \\min(r_t(\\theta) A_t, clip(r_t(\\theta), 1-\\epsilon, 1+\\epsilon) A_t) ]" - - -class KLPenaltyLoss(Loss): - r"""KL Penalty Loss - - KL Penalty to minimize: - - .. math:: L^{KL}(\theta) = \mathbb{E}[ KL( \pi_{\theta_{old}}(a_t | s_t) || \pi_{\theta}(a_t | s_t) ) ] - - where :math:`KL(.||.)` is the KL-divergence between two probability distributions. - """ - - def __init__(self): # p, q): - """ - Initialize the KL Penalty loss. - - Args: - p (torch.distributions.Distribution): 1st distribution - q (torch.distributions.Distribution): 2nd distribution - """ - super(KLPenaltyLoss, self).__init__() - # self.p = p - # self.q = q - - def compute(self, batch): - """ - Compute :math:`KL(p||q)`. - """ - curr_pi = batch.current['action_distributions'] - prev_pi = batch['action_distributions'] - - return torch.distributions.kl.kl_divergence(prev_pi, curr_pi) - - def latex(self): - return "\\mathbb{E}[ KL( \\pi_{\\theta_{old}}(a_t | s_t) || \\pi_{\\theta}(a_t | s_t) ) ]" - - -# class ForwardKLPenaltyLoss(KLPenaltyLoss): -# r"""Forward KL Penalty Loss""" -# pass -# -# -# class ReverseKLPenaltyLoss(KLPenaltyLoss): -# r"""Reverse KL Penalty Loss""" -# pass - - -class EntropyLoss(Loss): - r"""Entropy Loss - - Entropy loss, which is used to ensure sufficient exploration when maximized [1,2,3]: - - .. math:: L^{Entropy}(\theta) = H[ \pi_{\theta} ] - - where :math:`H[.]` is the Shannon entropy of the given probability distribution. - - References: - [1] "Simple Statistical Gradient-following Algorithms for Connectionist Reinforcement Learning", Williams, 1992 - [2] "Asynchronous Methods for Deep Reinforcement Learning", Mnih et al., 2016 - [3] "Proximal Policy Optimization Algorithms", Schulman et al., 2017 - """ - - def __init__(self): # approximator): - super(EntropyLoss, self).__init__() - - def compute(self, batch): - distribution = batch.current['action_distributions'] - entropy = distribution.entropy().mean() - return entropy - - # TODO: the functions defined here should be the same for `reward.py` and `loss.py`. Thus, we need to check numpy # TODO: or torch inside the functions. def ceil(x): @@ -763,10 +540,3 @@ def trunc(x): return y else: return torch.trunc(x) - - -# Tests -if __name__ == '__main__': - # compute the losses - loss = - FixedLoss(3) # + FixedLoss(2) - print(loss(2)) diff --git a/pyrobolearn/losses/losses.py b/pyrobolearn/losses/losses.py new file mode 100644 index 0000000..05d3b8c --- /dev/null +++ b/pyrobolearn/losses/losses.py @@ -0,0 +1,314 @@ +#!/usr/bin/env python +"""Defines the common loss functions that are used by the learning algorithm / optimizer. + +Losses are evaluated on model parameters, data batches / storages, or transitions tuples. +""" + +import torch + +from pyrobolearn.losses.loss import Loss + +__author__ = "Brian Delhaisse" +__copyright__ = "Copyright 2018, PyRoboLearn" +__credits__ = ["Brian Delhaisse"] +__license__ = "MIT" +__version__ = "1.0.0" +__maintainer__ = "Brian Delhaisse" +__email__ = "briandelhaisse@gmail.com" +__status__ = "Development" + + +class FixedLoss(Loss): + r"""Fixed Loss + + """ + def __init__(self, value): + super(FixedLoss, self).__init__() + self.value = value + + def compute(self, batch): + return self.value * batch + + +class L2Loss(Loss): + r"""L2 Loss + + Compute the L2 loss given by: :math:`1/2 * (y_{target} - y_{predict})^2` + """ + def __init__(self, target, approximator): + super(L2Loss, self).__init__() + self.target = target + self.approximator = approximator + + def compute(self, batch): + # based on approximator check what we need + return 0.5 * (self.target(batch) - self.approximator(batch)).pow(2).mean() + + +class ValueLoss(Loss): + r"""L2 loss for values + """ + def __init__(self): + super(ValueLoss, self).__init__() + + def compute(self, batch): + returns = batch['returns'] + values = batch.current['values'] + return 0.5 * (returns - values).pow(2).mean() + + +class HuberLoss(Loss): + r"""Huber Loss + + "In statistics, the Huber loss is a loss function used in robust regression, that is less sensitive to outliers + in data than the squared error loss." [1] + + This loss is given by [1]: + + .. math:: {\mathcal L}(\delta) = \left\{ \begin{array}{lc} 1/2 a^2 & for |a| \leq \delta, \\ + \delta (|a| - 1/2 \delta), & \mbox{otherwise} \left \end{array} + + "This function is quadratic for small values of :math:`a`, and linear for large values, with equal values and + slopes of the different sections at the two points where :math:`|a| = \delta`" [1]. + + In [2], this loss is used for DQN, where :math:`\delta=1`, and :math:`a` is the temporal difference error, that is, + :math:`a = Q(s,a) - (r + \gamma \max_a Q(s',a))` where :math:`(r + \gamma \max_a Q(s',a))` is the target function. + + References: + [1] Huber Loss (on Wikipedia): https://en.wikipedia.org/wiki/Huber_loss + [2] "Reinforcement Learning (DQN) Tutorial": + https://pytorch.org/tutorials/intermediate/reinforcement_q_learning.html + """ + def __init__(self, loss, delta=1.): + super(HuberLoss, self).__init__() + self.loss = loss + self.delta = delta + + def compute(self, batch): + a = self.loss(batch) + if abs(a) <= self.delta: + return 0.5 * torch.pow(a, 2) + return self.delta * (torch.abs(a) - 0.5 * self.delta) + + +class PGLoss(Loss): + r"""Policy Gradient Loss + + Compute the policy gradient loss which is maximized and given by: + + .. math:: L^{PG} = \mathbb{E}[ \log \pi_{\theta}(a_t | s_t) \psi_t ] + + where :math:`\psi_t` is the associated return estimator, which can be for instance, the total reward estimator + :math:`\psi_t = R(\tau)` (where :math:`\tau` represents the whole trajectory), the state action value estimator + :math:`\psi_t = Q(s_t, a_t)`, or the advantage estimator :math:`\psi_t = A_t = Q(s_t, a_t) - V(s_t)`. Other + estimators are also possible. + + The gradient with respect to the parameters :math:`\theta` is then given by: + + .. math:: g = \mathbb{E}[ \nabla_\theta \log \pi_{\theta}(a_t | s_t) + + References: + [1] "Proximal Policy Optimization Algorithms", Schulman et al., 2017 + [2] "High-Dimensional Continuous Control using Generalized Advantage Estimation", Schulman et al., 2016 + """ + + def __init__(self): + super(PGLoss, self).__init__() + + def compute(self, batch): + log_curr_pi = batch.current['action_distributions'] + log_curr_pi = log_curr_pi.log_probs(batch.current['actions']) + estimator = batch['estimator'] + loss = torch.exp(log_curr_pi) * estimator + return -loss.mean() + + def latex(self): + return "\\mathbb{E}[ r_t(\\theta) A_t ]" + + +class CPILoss(Loss): + r"""CPI Loss + + Conservative Policy Iteration objective which is maximized and defined in [1]: + + .. math:: L^{CPI}(\theta) = \mathbb{E}[ r_t(\theta) A_t ] + + where the expectation is taken over a finite batch of samples, :math:`A_t` is an estimator of the advantage fct at + time step :math:`t`, :math:`r_t(\theta)` is the probability ratio given by + :math:`r_t(\theta) = \frac{ \pi_{\theta}(a_t|s_t) }{ \pi_{\theta_{old}}(a_t|s_t) }`. + + References: + [1] "Approximately optimal approximate reinforcement learning", Kakade et al., 2002 + [2] "Proximal Policy Optimization Algorithms", Schulman et al., 2017 + """ + + def __init__(self): + """ + Initialize the CPI Loss. + """ + super(CPILoss, self).__init__() + + def compute(self, batch): # policy_distribution, old_policy_distribution, estimator): + # ratio = policy_distribution / old_policy_distribution + log_curr_pi = batch.current['action_distributions'] + log_curr_pi = log_curr_pi.log_probs(batch.current['actions']) + log_prev_pi = batch['action_distributions'] + log_prev_pi = log_prev_pi.log_probs(batch['actions']) + + ratio = torch.exp(log_curr_pi - log_prev_pi) + estimator = batch['estimator'] + + loss = ratio * estimator + return -loss.mean() + + def latex(self): + return "\\mathbb{E}[ r_t(\\theta) A_t ]" + + +class CLIPLoss(Loss): + r"""CLIP Loss + + Loss defined in [1] which is maximized and given by: + + .. math:: L^{CLIP}(\theta) = \mathbb{E}[ \min(r_t(\theta) A_t, clip(r_t(\theta), 1-\epsilon, 1+\epsilon) A_t) ] + + where the expectation is taken over a finite batch of samples, :math:`A_t` is an estimator of the advantage fct at + time step :math:`t`, :math:`r_t(\theta)` is the probability ratio given by + :math:`r_t(\theta) = \frac{ \pi_{\theta}(a_t|s_t) }{ \pi_{\theta_{old}}(a_t|s_t) }`. + + References: + [1] "Proximal Policy Optimization Algorithms", Schulman et al., 2017 + """ + + def __init__(self, clip=0.2): + """ + Initialize the loss. + + Args: + epsilon (float): clip parameter + """ + super(CLIPLoss, self).__init__() + self.eps = clip + + def compute(self, batch): # , policy_distribution, old_policy_distribution, estimator): + log_curr_pi = batch.current['action_distributions'] + log_curr_pi = log_curr_pi.log_probs(batch.current['actions']) + log_prev_pi = batch['action_distributions'] + log_prev_pi = log_prev_pi.log_probs(batch['actions']) + + ratio = torch.exp(log_curr_pi - log_prev_pi) + estimator = batch['estimator'] + + loss = torch.min(ratio * estimator, torch.clamp(ratio, 1.0-self.eps, 1.0+self.eps) * estimator) + return -loss.mean() + + def latex(self): + return "\\mathbb{E}[ \\min(r_t(\\theta) A_t, clip(r_t(\\theta), 1-\\epsilon, 1+\\epsilon) A_t) ]" + + +class KLPenaltyLoss(Loss): + r"""KL Penalty Loss + + KL Penalty to minimize: + + .. math:: L^{KL}(\theta) = \mathbb{E}[ KL( \pi_{\theta_{old}}(a_t | s_t) || \pi_{\theta}(a_t | s_t) ) ] + + where :math:`KL(.||.)` is the KL-divergence between two probability distributions. + """ + + def __init__(self): # p, q): + """ + Initialize the KL Penalty loss. + + Args: + p (torch.distributions.Distribution): 1st distribution + q (torch.distributions.Distribution): 2nd distribution + """ + super(KLPenaltyLoss, self).__init__() + # self.p = p + # self.q = q + + def compute(self, batch): + """ + Compute :math:`KL(p||q)`. + """ + curr_pi = batch.current['action_distributions'] + prev_pi = batch['action_distributions'] + + return torch.distributions.kl.kl_divergence(prev_pi, curr_pi) + + def latex(self): + return "\\mathbb{E}[ KL( \\pi_{\\theta_{old}}(a_t | s_t) || \\pi_{\\theta}(a_t | s_t) ) ]" + + +# class ForwardKLPenaltyLoss(KLPenaltyLoss): +# r"""Forward KL Penalty Loss""" +# pass +# +# +# class ReverseKLPenaltyLoss(KLPenaltyLoss): +# r"""Reverse KL Penalty Loss""" +# pass + + +class EntropyLoss(Loss): + r"""Entropy Loss + + Entropy loss, which is used to ensure sufficient exploration when maximized [1,2,3]: + + .. math:: L^{Entropy}(\theta) = H[ \pi_{\theta} ] + + where :math:`H[.]` is the Shannon entropy of the given probability distribution. + + References: + [1] "Simple Statistical Gradient-following Algorithms for Connectionist Reinforcement Learning", Williams, 1992 + [2] "Asynchronous Methods for Deep Reinforcement Learning", Mnih et al., 2016 + [3] "Proximal Policy Optimization Algorithms", Schulman et al., 2017 + """ + + def __init__(self): # approximator): + super(EntropyLoss, self).__init__() + + def compute(self, batch): + distribution = batch.current['action_distributions'] + entropy = distribution.entropy().mean() + return entropy + + +class MSBELoss(Loss): + r"""Mean-squared Bellman error + + The mean-squared Bellman error (MSBE) computes the Bellman error, also known as the one-step temporal difference + (TD). + + This is given by, in the case of the off-policy Q-learning TD algorithm: + .. math:: (r + \gamma (1-d) max_{a'} Q_{\phi}(s',a')) - Q_{\phi}(s,a), + + or in the case of the on-policy Sarsa TD algorithm: + .. math:: (r + \gamma (1-d) Q_{\phi}(s',a')) - Q_{\phi}(s,a) + + or in the case of the TD(0): + .. math:: (r + \gamma (1-d) V_{\phi}(s')) - V_{\phi}(s) + + where (r + \gamma (1-d) f(s,a)) is called (one-step return) the target. The target could also be, instead of the + one step return, the n-step return or the lambda-return value. [2] + These losses roughly tell us how closely the value functions satisfy the Bellman equation. Thus, by trying to + minimize them, we try to enforce the Bellman equations. + + References: + [1] https://spinningup.openai.com/en/latest/algorithms/ddpg.html + [2] "Reinforcement Learning: An Introduction", Sutton and Barto, 2018 + """ + + def __init__(self): + super(MSBELoss, self).__init__() + + def compute(self, batch): + pass + + +# Tests +if __name__ == '__main__': + # compute the losses + loss = - FixedLoss(3) # + FixedLoss(2) + print(loss(2)) diff --git a/pyrobolearn/optimizers/optimization.py b/pyrobolearn/optimizers/optimization.py new file mode 100644 index 0000000..e30f515 --- /dev/null +++ b/pyrobolearn/optimizers/optimization.py @@ -0,0 +1,75 @@ + + +class TrustRegion(object): + r"""Trust Region + + "Trust region is a term used in mathematical optimization to denote the subset of the region of the objective + function that is approximated using a model function (often a quadratic). If an adequate model of the objective + function is found within the trust region, then the region is expanded; conversely, if the approximation is poor, + then the region is contracted. Trust-region methods are also known as restricted-step methods. + + The fit is evaluated by comparing the ratio of expected improvement from the model approximation with the actual + improvement observed in the objective function. Simple thresholding of the ratio is used as the criterion for + expansion and contraction; a model function is "trusted" only in the region where it provides a reasonable + approximation. + + Trust-region methods are in some sense dual to line-search methods: trust-region methods first choose a step size + (the size of the trust region) and then a step direction, while line-search methods first choose a step direction + and then a step size." [1] + + References: + [1] https://en.wikipedia.org/wiki/Trust_region + """ + pass + + +class LineSearch(object): + r"""Line Search + + "In optimization, the line search strategy is one of two basic iterative approaches to find a local minimum + :math:`\mathbf{x}^*` of an objective function :math:`f:\mathbb{R}^{n} \to \mathbb{R}`. The other approach is trust + region. + + The line search approach first finds a descent direction along which the objective function :math:`f` will be + reduced and then computes a step size that determines how far :math:`\mathbf{x}` should move along that direction. + The descent direction can be computed by various methods, such as gradient descent, Newton's method and + Quasi-Newton method. The step size can be determined either exactly or inexactly. + + Here is an example gradient method that uses a line search in step 4. + + 1. Set iteration counter k = 0, and make an initial guess :math:`\mathbf{x}_{0}` for the minimum + 2. Repeat: + 3. Compute a descent direction :math:`\mathbf{p}_k` + 4. Choose :math:`\alpha_k` to 'loosely' minimize :math:`h(\alpha)=f(\mathbf{x}_k + \alpha \mathbf{p}_k)` over + :math:`\alpha \in \mathbb{R}_{+}` + 5. Update :math:`\mathbf{x}_{k+1} = \mathbf{x}_k + \alpha_k \mathbf{p}_k`, and :math:`k = k + 1` + 6. Until :math:`|| \nabla f( \mathbf{x}_k ) || < tolerance + + At the line search step (4) the algorithm might either exactly minimize :math:`h`, by solving + :math:`h'(\alpha _{k})=0`, or loosely, by asking for a sufficient decrease in :math:`h`. One example of the former + is conjugate gradient method. The latter is called inexact line search and may be performed in a number of ways, + such as a backtracking line search or using the Wolfe conditions. + + Like other optimization methods, line search may be combined with simulated annealing to allow it to jump over + some local minima." [1] + + References: + [1] https://en.wikipedia.org/wiki/Line_search + """ + pass + + +class BacktrackingLineSearch(LineSearch): + r"""Backtracking Line Search + + "In (unconstrained) minimization, a backtracking line search, a search scheme based on the Armijo-Goldstein + condition, is a line search method to determine the maximum amount to move along a given search direction. + It involves starting with a relatively large estimate of the step size for movement along the search direction, + and iteratively shrinking the step size (i.e., "backtracking") until a decrease of the objective function is + observed that adequately corresponds to the decrease that is expected, based on the local gradient of the + objective function." [1] + + References: + [1] https://en.wikipedia.org/wiki/Backtracking_line_search + """ + pass diff --git a/pyrobolearn/optimizers/torch_optimizer.py b/pyrobolearn/optimizers/torch_optimizer.py index 9900118..8687ed1 100644 --- a/pyrobolearn/optimizers/torch_optimizer.py +++ b/pyrobolearn/optimizers/torch_optimizer.py @@ -6,9 +6,11 @@ instantiating them, here we can pass the parameters at a later stage. References: [1] https://pytorch.org/docs/stable/optim.html + [2] https://github.com/sbarratt/torch_cg/tree/master """ # Pytorch optimizers +import torch import torch.nn as nn import torch.optim as optim @@ -200,3 +202,68 @@ class SGD(Optimizer): if self.max_grad_norm is not None: nn.utils.clip_grad_norm_(params, self.max_grad_norm) self.optimizer.step() + + +class CG(Optimizer): + r"""Conjugate Gradient + + "In mathematics, the conjugate gradient method is an algorithm for the numerical solution of particular systems + of linear equations, namely those whose matrix is symmetric and positive-definite. The conjugate gradient method + is often implemented as an iterative algorithm, applicable to sparse systems that are too large to be handled by + a direct implementation or other direct methods such as the Cholesky decomposition. Large sparse systems often + arise when numerically solving partial differential equations or optimization problems. + + Suppose we want to solve the system of linear equations :math:`\mathbf{A} \mathbf{x} = \mathbf{b}`, for the vector + :math:`\mathbf{x}`, where the known n x n matrix :math:`\mathbf{A}` is symmetric (i.e., + :math:`\mathbf{A}^\top = \mathbf{A}`), positive-definite (i.e. :math:`\mathbf{x}^\top \mathbf{Ax} > 0` for all + non-zero vectors :math:`x \in \mathbb{R}^n`), and real, and :math:`\mathbf{b}` is known as well. We denote the + unique solution of this system by :math:`\mathbf{x}^*`." [1] + + Notably, this can be used to solve :math:`\mathbf{Hx} = \mathbf{g}` where :math:`\mathbf{H}` is the Hessian + matrix, and :math:`g` is the gradient. + + References: + [1] Conjugate Gradient (Wikipedia): https://en.wikipedia.org/wiki/Conjugate_gradient_method + [2] Hessian matrix (Wikipedia): https://en.wikipedia.org/wiki/Hessian_matrix#Use_in_optimization + [3] Hessian-Vector products: https://justindomke.wordpress.com/2009/01/17/hessian-vector-products/ + [4] Torch CG: https://github.com/sbarratt/torch_cg + """ + + def __init__(self, threshold=1.e-8, max_iters=10, *args, **kwargs): + super(CG, self).__init__(*args, **kwargs) + self._threshold = threshold + self._max_iters = int(max_iters) + + def optimize(self, A, b, x=None): + """ + Return the solution :math:`x` to the system of linear equations :math:`Ax=b`. + + Notably, this can be used + + Args: + A (torch.Tensor): square real symmetric and positive-definite matrix. + b (torch.Tensor): known vector. + x (torch.Tensor, None): initial solution. If None, it will be the zero vector. + + Returns: + torch.Tensor: return the solution x to Ax=b. + """ + if x is None: + x = torch.zeros_like(b) + + r = b - A.matmul(x) + p = r + rr_old = r.t().matmul(r) + + for k in range(self._max_iters): + Ap = A.matmul(p) + alpha = rr_old / p.t().matmul(Ap) + x = x + alpha * p + r = r - alpha * Ap + rr_new = r.t().matmul(r) + if torch.sqrt(rr_new) < self._threshold: + break + p = r + rr_new / rr_old * p + rr_old = rr_new + + return x diff --git a/pyrobolearn/states/generators/state_generator.py b/pyrobolearn/states/generators/state_generator.py index 405550c..2240c4f 100644 --- a/pyrobolearn/states/generators/state_generator.py +++ b/pyrobolearn/states/generators/state_generator.py @@ -197,6 +197,23 @@ class QueueStateGenerator(StateGenerator): """Return the size of the queue.""" return self.qsize() + def __iter__(self): + """Return the iterator object itself.""" + self.cnt = 0 + return self + + def __next__(self): # only valid in Python 3 + """Return the next item in the sequence.""" + if self.cnt < self.qsize(): + self.cnt += 1 + return self.queue[self.cnt-1] + else: + raise StopIteration + + def next(self): # for Python 2 + """Return the next item in the sequence.""" + return self.__next__() + class FIFOQueueStateGenerator(QueueStateGenerator): r"""FIFO Queue Initial State Generator