# Patent

## Title
Directional optimization via EBW

## Abstract

An optimization system and method includes determining a best gradient as a sparse direction in a function having a plurality of parameters. The sparse direction includes a direction that maximizes change of the function. This maximum change of the function is determined by performing an optimization process that gives maximum growth subject to a sparsity regularized constraint. An extended Baum Welch (EBW) method can be used to identify the sparse direction. A best step size is determined along the sparse direction by finding magnitudes of entries of direction that maximizes the function restricted to the sparse direction. A solution is recursively refined for the function optimization using a processor and storage media.

## Background
BACKGROUND
1. Technical Field
The present invention relates to optimization of functions and more particularly to a system and method for optimizing functions with a large number of parameters.
2. Description of the Related Art
Optimization of functions with a large number of parameters is very complex. A few techniques for function optimization include an exact computation of a Hessian matrix for Newton type of optimization methods. This type of computation is difficult since it requires processing around a squared number of parameters. This results in a slow speed and memory overload in many instances.
Hill-climbing algorithms require step sizes being defined for each iteration and also have a relatively slow convergence. Optimization algorithms that are based on approximated versions of Hessians are counter-intuitive and difficult to customize to specific tasks in which a learning ratio is needed to be controlled.
Therefore, there is a need for an optimization method and system for optimization for a large number of parameters that is fast, intuitive and easy to customize to special tasks.
SUMMARY
An optimization system and method includes determining a best gradient as a sparse direction in a function having a plurality of parameters. The sparse direction includes a direction that maximizes change of the function. This maximum change of the function is determined by performing an optimization process that gives maximum growth subject to a sparsity regularized constraint. An extended Baum Welch (EBW) method can be used to identify the sparse direction. A best step size is determined along the sparse direction by finding magnitudes of entries of direction that maximizes the function restricted to the sparse direction. A solution is recursively refined for the function optimization using a processor and storage media.
An optimization method includes determining a best gradient as a sparse direction in a function having a plurality of parameters, the sparse direction including a direction that maximizes change of the function as determined by an optimization method where a sparse direction is defined from sparse regularized constraints. A best step size is determined along the sparse direction by finding magnitudes of entries of direction by solving an optimization problem that is constrained to the sparse direction. A solution to the optimization problem is recursively refined using a processor and storage media until a threshold constraint is met.
A system for optimization of a function with a large number of parameters includes a processor configured to determine a best gradient as a sparse direction in a function having a plurality of parameters. The sparse direction includes a direction that maximizes change of the function as determined by an optimization method where a sparse direction is defined from sparse regularized constraints. The processor is configured to determine a best step size along the sparse direction by solving an optimization problem to find magnitudes of entries of direction that are constrained to the sparse direction. A memory storage device includes a recursively refined stored solution for an optimization solution of the function.
These and other features and advantages will become apparent from the following detailed description of illustrative embodiments thereof, which is to be read in connection with the accompanying drawings.



BRIEF DESCRIPTION OF DRAWINGS
The disclosure will provide details in the following description of preferred embodiments with reference to the following figures wherein:
FIG. 1 is a block/flow diagram showing an optimization system and method in accordance with one illustrative embodiment;
FIG. 2 is a plot showing signal reconstruction with extended Extended Baum Welch (EBW) for 100,000 parameters;
FIG. 3 is a plot showing signal reconstruction with an Extended Kalman Filter (EKF);
FIG. 4 is a plot showing signal reconstruction with an Extended Baum Welch (EBW) method; and
FIG. 5 is a block/flow diagram showing an optimization system in accordance with the present principles.



DETAILED DESCRIPTION OF PREFERRED EMBODIMENTS
In accordance with the present principles, a system and method provides an optimization that can efficiently handle a large number of parameters and is simpler to use than conventional methods. In one embodiment, an optimization splits an optimization recursion process into parts, e.g., two parts. In a first part, optimizational direction is efficiently define via an Extended Baum-Welch (EBW) method for discrete parameters (EBWD). In a second part, a step along the direction is determined. A search for a best step along a direction can be represented as finding magnitudes of entries of direction by parallelizing entries into an independent set of parameters that can be processed independently of each other.
As will be appreciated by one skilled in the art, aspects of the present invention may be embodied as a system, method or computer program product. Accordingly, aspects of the present invention may take the form of an entirely hardware embodiment, an entirely software embodiment (including firmware, resident software, micro-code, etc.) or an embodiment combining software and hardware aspects that may all generally be referred to herein as a “circuit,” “module” or “system.” Furthermore, aspects of the present invention may take the form of a computer program product embodied in one or more computer readable medium(s) having computer readable program code embodied thereon.
Any combination of one or more computer readable medium(s) may be utilized. The computer readable medium may be a computer readable signal medium or a computer readable storage medium. A computer readable storage medium may be, for example, but not limited to, an electronic, magnetic, optical, electromagnetic, infrared, or semiconductor system, apparatus, or device, or any suitable combination of the foregoing. More specific examples (a non-exhaustive list) of the computer readable storage medium would include the following: an electrical connection having one or more wires, a portable computer diskette, a hard disk, a random access memory (RAM), a read-only memory (ROM), an erasable programmable read-only memory (EPROM or Flash memory), an optical fiber, a portable compact disc read-only memory (CD-ROM), an optical storage device, a magnetic storage device, or any suitable combination of the foregoing. In the context of this document, a computer readable storage medium may be any tangible medium that can contain, or store a program for use by or in connection with an instruction execution system, apparatus, or device.
A computer readable signal medium may include a propagated data signal with computer readable program code embodied therein, for example, in baseband or as part of a carrier wave. Such a propagated signal may take any of a variety of forms, including, but not limited to, electro-magnetic, optical, or any suitable combination thereof. A computer readable signal medium may be any computer readable medium that is not a computer readable storage medium and that can communicate, propagate, or transport a program for use by or in connection with an instruction execution system, apparatus, or device.
Program code embodied on a computer readable medium may be transmitted using any appropriate medium, including but not limited to wireless, wireline, optical fiber cable, RF, etc., or any suitable combination of the foregoing. Computer program code for carrying out operations for aspects of the present invention may be written in any combination of one or more programming languages, including an object oriented programming language such as Java, Smalltalk, C++ or the like and conventional procedural programming languages, such as the “C” programming language or similar programming languages. The program code may execute entirely on the user's computer, partly on the user's computer, as a stand-alone software package, partly on the user's computer and partly on a remote computer or entirely on the remote computer or server. In the latter scenario, the remote computer may be connected to the user's computer through any type of network, including a local area network (LAN) or a wide area network (WAN), or the connection may be made to an external computer (for example, through the Internet using an Internet Service Provider).
Aspects of the present invention are described below with reference to flowchart illustrations and/or block diagrams of methods, apparatus (systems) and computer program products according to embodiments of the invention. It will be understood that each block of the flowchart illustrations and/or block diagrams, and combinations of blocks in the flowchart illustrations and/or block diagrams, can be implemented by computer program instructions. These computer program instructions may be provided to a processor of a general purpose computer, special purpose computer, or other programmable data processing apparatus to produce a machine, such that the instructions, which execute via the processor of the computer or other programmable data processing apparatus, create means for implementing the functions/acts specified in the flowchart and/or block diagram block or blocks.
These computer program instructions may also be stored in a computer readable medium that can direct a computer, other programmable data processing apparatus, or other devices to function in a particular manner, such that the instructions stored in the computer readable medium produce an article of manufacture including instructions which implement the function/act specified in the flowchart and/or block diagram block or blocks.
The computer program instructions may also be loaded onto a computer, other programmable data processing apparatus, or other devices to cause a series of operational steps to be performed on the computer, other programmable apparatus or other devices to produce a computer implemented process such that the instructions which execute on the computer or other programmable apparatus provide processes for implementing the functions/acts specified in the flowchart and/or block diagram block or blocks.
The flowchart and block diagrams in the FIGS. illustrate the architecture, functionality, and operation of possible implementations of systems, methods and computer program products according to various embodiments of the present invention. In this regard, each block in the flowchart or block diagrams may represent a module, segment, or portion of code, which comprises one or more executable instructions for implementing the specified logical function(s). It should also be noted that, in some alternative implementations, the functions noted in the block may occur out of the order noted in the figures. For example, two blocks shown in succession may, in fact, be executed substantially concurrently, or the blocks may sometimes be executed in the reverse order, depending upon the functionality involved. It will also be noted that each block of the block diagrams and/or flowchart illustration, and combinations of blocks in the block diagrams and/or flowchart illustration, can be implemented by special purpose hardware-based systems that perform the specified functions or acts, or combinations of special purpose hardware and computer instructions.
Referring now to the drawings in which like numerals represent the same or similar elements and initially to FIG. 1, a block/flow diagram illustratively shows a system/method for optimizing a function in accordance with one exemplary embodiment. In block 10, a function to be optimized is provided. The function may include a large number of parameters (e.g., hundreds or even thousands) which may provide one or more directions which make the solution very complex. A direction of the function at a point is defined as a weighted sum of base vectors that belongs to a linear subspace that passes through this point.
In block 12, an optimization recursion process is split into at least two parts. In block 14, optimizational direction (sparse direction) is efficiently defined via an optimization method, e.g., an Extended Baum-Welch (EBW) method for discrete parameters (EBWD). This will be explained in greater detail below. In block 16, a step along the sparse direction is determined. A search for a best step along directions can be represented as finding magnitudes of entries that are weights identified as non-zero components in the sparse direction. In block 18, an optimized parameter set or optimized data are provided or output for the given function.
It should be understood that the present embodiments may be implemented in a plurality of technologies for a plurality of different optimization applications. As illustratively described the present principles are described with reference to a general function; however, the function may represent any of a plurality of quantities. For example, the function may represent a scheduling function with a plurality of parametric constraints, may represent a vector or tensor with a large number of variables, may represent the motion of an object three dimensions, etc. Other examples of the function include, e.g., a maximum mutual information function that represents information between data and labels of this data, a maximum entropy function, functions associated with a support vector machine, etc. This list is non-exhaustive as it represents but a few possible applications which can be described by a function to be optimized.
In block 14, a direction is defined. EBWD was initially suggested for optimizing an arbitrary (differential) function over a probability discrete domain. To apply EBWD to an arbitrary function for unconstrained optimization, one can represent directions as the following. Let F(x) be a differentiable function (block 10 of FIG. 1) of




  
    
      
        
          
            x
            =
            
              {
              
                x
                i
              
              }
            
          
          ,
          
            v
            =
            
              
                {
                
                  v
                  i
                
                }
              
              ∈
              
                R
                n
              
            
          
          ,
          
            i
            =
            1
          
          ,
          
            ...
            ⁢
            
              
            
            ⁢
            n
          
          ,
          where
        
        ⁢
        
          

        
        ⁢
        
          
            
              max
              v
            
            ⁢
            
              F
              ⁡
              
                (
                
                  
                    v
                    0
                  
                  +
                  v
                
                )
              
            
          
          -
          
            F
            ⁡
            
              (
              
                v
                0
              
              )
            
          
        
        ⁢
        
          

        
        ⁢
        
          
            subject
            ⁢
            
              
            
            ⁢
            to
            ⁢
            
              
            
            ⁢
            
              (
              
                s
                .
                t
                .
              
              )
            
            ⁢
            
              
                
                v
                
              
              1
            
          
          <
          ɛ
        
      
    
    
      
        (
        1
        )
      
    
  




R is one dimensional real space; Rn is n dimensional real space. Direction is a vector v={vi} consisting of variables vi. ∥v∥1 is norm 1 of a vector v=(v1, v2, . . . vL) and is defined as





  
    
      
      v
      
    
    1
  
  =
  
    
      ∑
      
        i
        =
        1
      
      L
    
    ⁢
    
      
      
        v
        i
      
      
    
  




where vi are real numbers (components of the vector v) and |vi| are absolute values of vi. The direction vector v belongs to a subspace that passes through the point v0.

For sufficiently small ε, this approximates finding a best gradient:




  
    
      
        
          max
          v
        
        ⁢
        
          
            
              ∂
              
                F
                ⁡
                
                  (
                  
                    
                      v
                      0
                    
                    +
                    v
                  
                  )
                
              
            
            
              ∂
              v
            
          
          .
        
      
    
    
      
        (
        2
        )
      
    
  




Therefore, the gradient of function F(v) at point v0 is computed as





  
    
      
        
          d
          F
        
        
          d
          v
        
      
      
    
    
      v
      0
    
  
  ⁢
  
    
      (
      v
      )
    
    .
  




We can replace this expression with the following approximation:





  
    
      F
      ⁡
      
        (
        
          
            v
            0
          
          +
          v
        
        )
      
    
    -
    
      F
      ⁡
      
        (
        
          v
          +
          w
        
        )
      
    
  
  v




adding the constraint that ∥v∥1<ε. Essentially we replace the step of finding a gradient of maximum growth (the usual step in optimization) with a problem of finding a gradient that gives maximum growth subject to a sparsity regularized constraint. So the optimization is done in accordance with this sparse constraint, e.g., ∥v∥1<ε. Thus, the gradient problem becomes Eq. (1).

A ball for a norm 1 is a set of vectors that satisfies inequality ∥v∥1<ε where ε is some number (e.g., a threshold). The smaller epsilon, the smaller the ball is. The ball is a metaphoric name. For example, a best gradient includes identifying a size ε of a directional ball for determining the best direction for a directional optimization problem in parameter space of the plurality of parameters. Given the sparse direction, a larger directional ball R along the sparse direction may be determined (e.g, by performing the extended Baum-Welch method) so that a sub-optimal size inside R can be identified. The R ball is just a notation of a different ball.
In block 16, magnitudes are defined. After a best direction {tilde over (v)} is defined from Eq. (1), we need to find xεRn that optimizes maxs F(v0+s{tilde over (v)}) (3) where sv={si}{vi}={sivi} is a multiplication by components.
This can be performed by either using traditional optimization methods or formulating this as a constrained problem with linear domain constraints that can be solved via EBWD. EBWD is a method for optimizing a general function of a discrete probability domain. The method uses a first derivative and projection on the probability domain. Therefore, it is a large scale method.
Then, a new point can be defined x=v0+sv as a result of the phases performed in blocks 14 and 16. This is an inteimmediate result. To get a final result blocks 14 and 16 need to be performed recursively, updating point x until some threshold condition(s) is met.
In block 12, the recursion includes starting with a point xk, k=0. A sparse direction {tilde over (v)} is computed using Eq. (1) (block 14). A best step s={si} is computed using Eq. (3), and then next point xk=s*{tilde over (v)} (block 16). If threshold criteria are met, stop the process and consider xk as the outcome of the optimization process. Otherwise, k=k+1 and go to block 14.
Using some mathematical manipulations, we can rewrite function F(v) into another space G(z) and thus the above gradient computation is formulated as



  
    
      max
      z
    
    ⁢
    
      G
      ⁡
      
        (
        z
        )
      
    
  




  
    
      s
      .
      t
      .
      
        
      
      ⁢
      
        
          
          z
          
        
        q
      
    
    ⁢
    ɛ
  



Norm q, where q is some positive number, of a vector z=(z1, z2, . . . zL) is defined as




  
    
      
      z
      
    
    q
  
  =
  
    
      ∑
      
        i
        =
        1
      
      L
    
    ⁢
    
      
        
          (
          
            
              
              
                z
                i
              
              
            
            q
          
          )
        
        
          1
          /
          q
        
      
      .
    
  




This problem cannot be solved using available algorithms, such as, e.g., LASSO. One problem with LASSO is that it does not work for a large number of parameters. Therefore, we can take F(v) and again with some mathematical manipulation, express it as a function G(z) and thus the gradient computation is formulated as:




  
    
      max
      z
    
    ⁢
    
      G
      ⁡
      
        (
        z
        )
      
    
  




  
    
      
        s
        .
        t
        .
        
          
        
        ⁢
        
          ∑
          
            z
            i
          
        
      
      =
      1
    
    ,
    
      
        z
        1
      
      ≥
      0
    
  



This solution has a simple EBW solution for finding magnitudes employed in determined the step sized along a direction. The computation of the gradient can be represented using EBW. EBW for finding magnitude may include the following. Let ε>0 be an upper boundary that ∥v∥1 does not exceed. Let us re-formulate (3) as the following




  
    
      
        
          
            max
            m
          
          ⁢
          
            F
            ⁡
            
              (
              
                
                  x
                  0
                
                +
                xv
              
              )
            
          
        
        ⁢
        
          

        
        ⁢
        
          
            s
            .
            t
            .
            
              
            
            ⁢
            
              
                
                x
                
              
              1
            
          
          ≤
          ɛ
        
      
    
    
      
        (
        4
        )
      
    
  




The solution to Eq. (4) is given below.

DISCRETE EBW: Let F(x) be a differentiable function of x={xi}εRn, i=1, . . . n . Let us consider the following problem:




  
    
      
        
          max
          ⁢
          
            
          
          ⁢
          
            F
            ⁡
            
              (
              x
              )
            
          
        
        ⁢
        
          

        
        ⁢
        
          
            s
            .
            t
            .
            
              
            
            ⁢
            
              
                
                x
                
              
              1
            
          
          ≤
          ɛ
        
      
    
    
      
        (
        5
        )
      
    
  




We solve this problem by transforming Eq. (26), as set forth below, into a problem over a probability domain for which EBW update rules exist. Let us transform this problem into the following: Turning (26) into an equality, let x0≧0 be a dummy variable and let us consider the following problem with the dummy variable added:




  
    
      
        
          max
          ⁢
          
            
          
          ⁢
          
            F
            ⁡
            
              (
              x
              )
            
          
        
        ⁢
        
          

        
        ⁢
        
          
            
              subject
              ⁢
              
                
              
              ⁢
              
                
                  
                  x
                  
                
                1
              
            
            +
            
              x
              0
            
          
          ≤
          
            ɛ
            .
          
        
      
    
    
      
        (
        6
        )
      
    
  




Note adding a dumb non-negative variable is employed to map a fractional ball constrained inequality into a fractionally constrained multi-dimensional sphere. Turning ε into 1 (norm 1), let vi=xi/ε, i=0, . . . n . Let F(x)=F({εvi})=G(v) let us consider the following problem





  
    
      
        
          max
          ⁢
          
            
          
          ⁢
          
            G
            ⁡
            
              (
              v
              )
            
          
        
        ⁢
        
          

        
        ⁢
        
          
            
              
                subject to
              
              ⁢
              
                
              
              ⁢
              
                
                  ∑
                  i
                
                ⁢
                
                  
                    
                    
                      v
                      i
                    
                    
                  
                  1
                
              
            
            +
            
              v
              0
            
          
          =
          1.
        
      
    
    
      
        (
        7
        )
      
    
  




This transforms the fractional multi-dimensional sphere into a norm one sphere. The norm one sphere can be mapped into a simplex of non-negative numbers, and an objective function over the norm one sphere may be transformed into the objective function over the simplex of non-negative numbers.
Reduction to discrete EBW: Let us consider the following problem:




  
    
      
        
          max
          ⁢
          
            
          
          ⁢
          
            F
            ⁡
            
              (
              x
              )
            
          
        
        ⁢
        
          

        
        ⁢
        
          
            subject
            ⁢
            
              
            
            ⁢
            to
            ⁢
            
              
            
            ⁢
            
              
                
                x
                
              
              1
            
          
          =
          1
        
      
    
    
      
        (
        8
        )
      
    
  




Let us set




  
    
      
        
          
            
            x
            
          
          1
        
        =
        
          
            ∑
            
              
              
                x
                i
              
              
            
          
          =
          
            
              ∑
              i
            
            ⁢
            
              
                σ
                ⁡
                
                  (
                  
                    x
                    i
                  
                  )
                
              
              ⁢
              
                x
                i
              
            
          
        
      
    
    
      
        (
        9
        )
      
    
  




where σ(xi)=1 if xi≧0 and σ(xi)=−1 if xi≦0. Let

G(z)=F({σ(xi)σ(xi)xi})=F(σ(xi)zi)  (10)

where z={zi}={σ(xi)xi} Then the Eq. (26) is equivalent to the following problem:





  
    
      
        
          max
          ⁢
          
            
          
          ⁢
          
            G
            ⁡
            
              (
              z
              )
            
          
        
        ⁢
        
          

        
        ⁢
        
          
            
              subject
              ⁢
              
                
              
              ⁢
              to
              ⁢
              
                
              
              ⁢
              
                ∑
                
                  z
                  i
                
              
            
            =
            1
          
          ,
          
            
              z
              i
            
            ≥
            0.
          
        
      
    
    
      
        (
        11
        )
      
    
  




This maps the sphere normalized to one into a simplex of non-negative numbers. An objective function over the norm one sphere may be transformed into an objective function over the simplex.
Problems of the type of Equation (11) have EBW based solutions that can be described as the following:




  
    
      
        
          z
          i
          
            t
            +
            1
          
        
        =
        
          
            
              z
              i
              
                t
                +
                1
              
            
            ⁡
            
              (
              D
              )
            
          
          =
          
            
              
                (
                
                  
                    c
                    i
                  
                  +
                  D
                
                )
              
              ⁢
              
                z
                i
                t
              
            
            
              
                
                  ∑
                  j
                
                ⁢
                
                  
                    c
                    j
                  
                  ⁢
                  
                    z
                    j
                    t
                  
                
              
              +
              D
            
          
        
      
    
    
      
        (
        12
        )
      
    
  




where t is an index and





  
    
      
        
          c
          j
        
        =
        
          
            
              c
              j
            
            ⁡
            
              (
              
                z
                t
              
              )
            
          
          =
          
            
              
                ∂
                
                  G
                  ⁡
                  
                    (
                    
                      z
                      t
                    
                    )
                  
                
              
              
                ∂
                
                  z
                  j
                  t
                
              
            
            .
          
        
      
    
    
      
        (
        13
        )
      
    
  




It has been proved that Eq. (29) provides growth transformations for sufficiently large D if G is a differentiable function in z, i.e., consequent applications of Eq. (29) converge to a local maximum of G(z) in the domain Σzi=1, zi≧0. D is a constant that controls growth of the function that is being optimized.
Reduction to discrete EBW is illustratively depicted in a recursive form. Eq. (8) is represented below in a recursive form.

I) set initial conditions for the recursion. Initially, x0={xi0εR1}s.t.∥x0∥1=1  (14).

II) zit=σ(xit)xit  (15).

III) Gt(zt)=F({σ(xit)σ(xit)xit})=F(σ(xit)zit)  (16)

(substitution from block 24).




  
    
      
        
          
            IV
            )
          
          ⁢
          
            
          
          ⁢
          
            c
            j
          
        
        =
        
          
            
              c
              j
            
            ⁡
            
              (
              
                z
                t
              
              )
            
          
          =
          
            
              
                ∂
                
                  
                    G
                    t
                  
                  ⁡
                  
                    (
                    
                      z
                      t
                    
                    )
                  
                
              
              
                ∂
                
                  z
                  j
                  t
                
              
            
            .
          
        
      
    
    
      
        (
        17
        )
      
    
  
  
    
      
        
          
            V
            )
          
          ⁢
          
            
          
          ⁢
          
            z
            i
            
              t
              +
              1
            
          
        
        =
        
          
            
              z
              i
              
                t
                +
                1
              
            
            ⁡
            
              (
              D
              )
            
          
          =
          
            
              
                
                  (
                  
                    
                      c
                      i
                    
                    +
                    D
                  
                  )
                
                ⁢
                
                  z
                  i
                  t
                
              
              
                
                  
                    ∑
                    j
                  
                  ⁢
                  
                    
                      c
                      j
                    
                    ⁢
                    
                      z
                      j
                      t
                    
                  
                
                +
                D
              
            
            .
          
        
      
    
    
      
        (
        18
        )
      
    
  
  
    
      
        
          
            VI
            )
          
          ⁢
          
            
          
          ⁢
          
            D
            *
          
        
        =
        
          arg
          ⁢
          
            
          
          ⁢
          
            
              max
              D
            
            ⁢
            
              
                
                  G
                  t
                
                ⁡
                
                  (
                  
                    {
                    
                      
                        z
                        i
                        
                          t
                          +
                          1
                        
                      
                      ⁡
                      
                        (
                        D
                        )
                      
                    
                    }
                  
                  )
                
              
              .
            
          
        
      
    
    
      
        (
        19
        )
      
    
  
  
    
      
        
          
            VII
            )
          
          ⁢
          
            
          
          ⁢
          
            x
            i
            
              t
              +
              1
            
          
        
        =
        
          
            
              z
              i
              
                t
                +
                1
              
            
            ⁡
            
              (
              
                D
                *
              
              )
            
          
          .
        
      
    
    
      
        (
        20
        )
      
    
  




VIII) The program returns to step II. The quantities computed are indexed using index i.
Special numerical example: Let us start with a simple simulation example of a cubic homogenous polynomial




  
    
      
        P
        =
        
          
            P
            ⁡
            
              (
              z
              )
            
          
          =
          
            
              ∑
              υ
            
            ⁢
            
              
                a
                υ
              
              ⁢
              
                z
                υ
              
            
          
        
      
    
    
      
        (
        21
        )
      
    
  




where υ={υ0, υ1, υ2|υi≧0, Συi=3} is a tuple of integers and





  
    z
    υ
  
  =
  
    
      ∏
      
        i
        =
        0
      
      
        i
        =
        2
      
    
    ⁢
    
      
        z
        i
        
          υ
          i
        
      
      .
    
  




In a case that Gt(zt)=P(zt):





  
    
      
        
          
            ∂
            
              
                G
                t
              
              ⁡
              
                (
                
                  z
                  t
                
                )
              
            
          
          
            ∂
            
              z
              j
              t
            
          
        
        =
        
          ∑
          
            
              a
              υ
            
            ⁢
            
              υ
              j
            
            ⁢
            
              z
              
                
                  υ
                  _
                
                ⁡
                
                  (
                  j
                  )
                
              
            
          
        
      
    
    
      
        (
        22
        )
      
    
  




where υ(j)i=υi if i≠j and υ(j)j=υj−1.

CONTINUOUS CASE (GAUSSIAN): Let us start with a simple simulation example of cubic homogenous polynomials

P=P(Z)=Σaυzυ  (23)

Q=Q(Z)=bυzυ  (24)

with non-negative coefficients aυ, bυ where

zi(x, μ, σ)=1/2 π1/2 exp(−1/2 (xi−μ)2/σ2  (25)

iε{0,1,2}, xi are arbitrary scalars (we do not require them to be nonnegative or sum up to 1) and μ, σ are model parameters. υ={υ0, υ1, υ2|υi≧0, Συi=3} is a tuple of integers and





  
    z
    υ
  
  =
  
    
      ∏
      
        i
        =
        0
      
      
        i
        =
        2
      
    
    ⁢
    
      
        z
        i
        
          υ
          i
        
      
      .
    
  




Next, let F(z)=P(z)/Q(z) and z={zi}.

Let us consider the following problem:




  
    
      
        
          max
          ⁢
          
            
          
          ⁢
          
            F
            ⁡
            
              (
              x
              )
            
          
        
        ⁢
        
          

        
        ⁢
        
          
            subject
            ⁢
            
              
            
            ⁢
            
              
                
                x
                
              
              1
            
          
          =
          1
        
      
    
    
      
        (
        26
        )
      
    
  




We can suggest two different solutions for this problem. The description of the solution is as follows. Let t=0, 1, . . . .

G(zt)=F({σ(xit)zit)  (27)

where zi0=σ(xi)xi, zit=σ(xit)xit and xit, zit are defined recursively as follows:





  
    
      
        
          c
          j
          t
        
        =
        
          
            
              c
              j
              t
            
            ⁡
            
              (
              
                z
                t
              
              )
            
          
          =
          
            
              
                ∂
                log
              
              ⁢
              
                
              
              ⁢
              
                G
                ⁡
                
                  (
                  
                    z
                    t
                  
                  )
                
              
            
            
              ∂
              
                z
                j
                t
              
            
          
        
      
    
    
      
        (
        28
        )
      
    
  
  
    
      
        
          x
          i
          
            t
            +
            1
          
        
        =
        
          
            
              x
              i
              
                t
                +
                1
              
            
            ⁡
            
              (
              
                D
                ^
              
              )
            
          
          =
          
            
              
                
                  (
                  
                    
                      c
                      i
                      t
                    
                    +
                    
                      D
                      ^
                    
                  
                  )
                
                ⁢
                
                  z
                  i
                  t
                
              
              
                
                  
                    ∑
                    j
                  
                  ⁢
                  
                    
                      c
                      j
                      t
                    
                    ⁢
                    
                      z
                      j
                      t
                    
                  
                
                +
                
                  D
                  ^
                
              
            
            .
          
        
      
    
    
      
        (
        29
        )
      
    
  




Here {circumflex over (D)} can be chosen as




  
    
      
        
          D
          ^
        
        =
        
          arg
          ⁢
          
            
          
          ⁢
          
            
              max
              D
            
            ⁢
            
              
                F
                ⁡
                
                  (
                  
                    
                      x
                      i
                      
                        t
                        +
                        1
                      
                    
                    ⁡
                    
                      (
                      D
                      )
                    
                  
                  )
                
              
              .
            
          
        
      
    
    
      
        (
        30
        )
      
    
  




D is a constant that controls growth of the function that is being optimized when {circumflex over (D)} is updated. {circumflex over (D)} is chosen as which D gives maximum growth (e.g., this needs a solution of an optimization problem (Eq. 30) in D to find {circumflex over (D)}).
In what follows, we give more detailed calculations for Eq. (28). Let




  
    
      
        z
        =
        
          
            
              z
              i
              t
            
            ⁡
            
              (
              
                x
                ,
                μ
                ,
                σ
              
              )
            
          
          =
          
            
              1
              /
              2
            
            ⁢
            
              π
              
                1
                /
                2
              
            
            ⁢
            
              exp
              ⁡
              
                (
                
                  
                    
                      -
                      1
                    
                    /
                    2
                  
                  ⁢
                  
                    
                      
                        (
                        
                          
                            x
                            i
                            t
                          
                          -
                          μ
                        
                        )
                      
                      2
                    
                    /
                    
                      σ
                      2
                    
                  
                
                )
              
            
          
        
      
    
    
      
        (
        31
        )
      
    
  
  
    
      
        
          
            Let
            ⁢
            
              
            
            ⁢
            
              H
              ⁡
              
                (
                
                  {
                  
                    
                      z
                      i
                      t
                    
                    ⁡
                    
                      (
                      
                        x
                        ,
                        μ
                        ,
                        σ
                      
                      )
                    
                  
                  }
                
                )
              
            
          
          =
          
            G
            ⁡
            
              (
              
                {
                
                  x
                  i
                  t
                
                }
              
              )
            
          
        
        ⁢
        
          

        
        ⁢
        Then
        ,
        
          
            
              
                ∂
                log
              
              ⁢
              
                
              
              ⁢
              
                G
                ⁡
                
                  (
                  
                    x
                    t
                  
                  )
                
              
            
            
              ∂
              
                x
                i
                t
              
            
          
          =
          
            
              
                
                  
                    ∂
                    log
                  
                  ⁢
                  
                    
                  
                  ⁢
                  
                    H
                    ⁡
                    
                      (
                      
                        z
                        t
                      
                      )
                    
                  
                
                
                  ∂
                  
                    z
                    i
                    t
                  
                
              
              ⁢
              
                
                  ∂
                  
                    z
                    i
                    t
                  
                
                
                  ∂
                  
                    x
                    i
                    t
                  
                
              
            
            =
            
              G
              ⁡
              
                (
                
                  {
                  
                    x
                    i
                    t
                  
                  }
                
                )
              
            
          
        
        ,
        where
      
    
    
      
        (
        32
        )
      
    
  
  
    
      
        
          h
          i
          t
        
        =
        
          
            
              h
              i
              t
            
            ⁡
            
              (
              z
              )
            
          
          =
          
            
              
                
                  ∂
                  log
                
                ⁢
                
                  
                
                ⁢
                
                  R
                  ⁡
                  
                    (
                    z
                    )
                  
                
              
              
                ∂
                
                  z
                  i
                
              
            
            =
            
              
                
                  
                    ∂
                    
                      P
                      ⁡
                      
                        (
                        z
                        )
                      
                    
                  
                  
                    ∂
                    
                      z
                      i
                    
                  
                
                /
                
                  P
                  ⁡
                  
                    (
                    z
                    )
                  
                
              
              =
              
                
                  
                    ∂
                    
                      Q
                      ⁡
                      
                        (
                        z
                        )
                      
                    
                  
                  
                    ∂
                    
                      z
                      i
                    
                  
                
                /
                
                  
                    Q
                    ⁡
                    
                      (
                      z
                      )
                    
                  
                  .
                
              
            
          
        
      
    
    
      
        (
        33
        )
      
    
  




This implies:

cit=xi└hitzit(−xit+μ)/σ2+C┘  (34)

In simulation experiments one can choose xi randomly and derive cit from Eq. (34) and update yit by solving Eq. (29). C is a tuning parameter that guarantees convergence. We need to check that updates in Eq. (34) and (29) converge to a local maximum of G(x) for sufficiently large C. Specifically, check that the following inequality is fulfilled for sufficiently large C: G(xt+1)≧G(xt). One can also find an optimal C in simulation experiments.

We consider the following unconstrained problem:

maxQ(P, ∥P∥1)  (35)

where PεRL and RL is L dimensional real space.

General loop: Eq. (35) can be solved via the following steps:

    
    
        1. Set initial value P0 
        2. Set d0=∥P0∥1 
        3. Find Pt+1=Pt+1(dt) in a functional form of dt s.t. ∥Pt+1∥1=dt and

Q(Pt+1,dt)≧Q(Pt,dt)  (36)

        4. Find dt+1 s.t.

Q(Pt+1(dt+1),dt+1)≧Q(Pt+1,dt)  (37)

        5. Go to Step 3.
    
    


The solution of Step 3 includes: Set x0=Pt+1/dt and F(x0)=Q(Pt+1,dt) Run several “EBW” recursions:




  
    
      
        
          
            a
            )
          
          ⁢
          
            
          
          ⁢
          Set
          ⁢
          
            
          
          ⁢
          initial
          ⁢
          
            
          
          ⁢
          conditions
          ⁢
          
            :
          
          ⁢
          
            
          
          ⁢
          
            x
            0
          
        
        =
        
          
            
              {
              
                
                  x
                  i
                  0
                
                ∈
                
                  R
                  1
                
              
              }
            
            ⁢
            
              s
              .
              t
              .
              
                
                  
                  
                    x
                    0
                  
                  
                
                1
              
            
          
          =
          1
        
      
    
    
      
        (
        38
        )
      
    
  
  
    
      
        
          
            b
            )
          
          ⁢
          
            
          
          ⁢
          
            z
            i
            t
          
        
        =
        
          
            σ
            ⁡
            
              (
              
                x
                i
                t
              
              )
            
          
          ⁢
          
            x
            i
            t
          
        
      
    
    
      
        (
        39
        )
      
    
  
  
    
      
        
          
            c
            )
          
          ⁢
          
            
          
          ⁢
          
            
              G
              t
            
            ⁡
            
              (
              
                z
                t
              
              )
            
          
        
        =
        
          
            F
            ⁡
            
              (
              
                {
                
                  
                    σ
                    ⁡
                    
                      (
                      
                        x
                        i
                        t
                      
                      )
                    
                  
                  ⁢
                  
                    σ
                    ⁡
                    
                      (
                      
                        x
                        i
                        t
                      
                      )
                    
                  
                  ⁢
                  
                    x
                    i
                    t
                  
                
                }
              
              )
            
          
          =
          
            F
            ⁡
            
              (
              
                
                  σ
                  ⁡
                  
                    (
                    
                      x
                      i
                      t
                    
                    )
                  
                
                ⁢
                
                  z
                  i
                  t
                
              
              )
            
          
        
      
    
    
      
        (
        40
        )
      
    
  
  
    
      
        
          
            d
            )
          
          ⁢
          
            
          
          ⁢
          
            c
            j
          
        
        =
        
          
            
              c
              j
            
            ⁡
            
              (
              
                z
                t
              
              )
            
          
          =
          
            
              ∂
              
                
                  G
                  t
                
                ⁡
                
                  (
                  
                    z
                    t
                  
                  )
                
              
            
            
              ∂
              
                z
                j
                t
              
            
          
        
      
    
    
      
        (
        41
        )
      
    
  
  
    
      
        
          
            e
            )
          
          ⁢
          
            
          
          ⁢
          
            z
            i
            
              t
              +
              1
            
          
        
        =
        
          
            
              z
              i
              
                t
                +
                1
              
            
            ⁡
            
              (
              D
              )
            
          
          =
          
            
              
                (
                
                  
                    c
                    i
                  
                  +
                  D
                
                )
              
              ⁢
              
                z
                i
                t
              
            
            
              
                
                  ∑
                  j
                
                ⁢
                
                  
                    c
                    j
                  
                  ⁢
                  
                    z
                    j
                    t
                  
                
              
              +
              D
            
          
        
      
    
    
      
        (
        42
        )
      
    
  




For some applications we can also use an approximation form of Eq. (29).





  
    
      
        
          z
          i
          
            t
            +
            1
          
        
        =
        
          
            
              z
              i
              
                t
                +
                1
              
            
            ⁡
            
              (
              D
              )
            
          
          =
          
            
              z
              i
              t
            
            +
            
              
                
                  ∑
                  
                    
                      c
                      i
                    
                    ⁢
                    
                      z
                      i
                      t
                    
                  
                
                D
              
              ⁢
              
                
                  (
                  
                    1
                    -
                    
                      
                        
                          z
                          i
                          t
                        
                        ⁢
                        
                          c
                          i
                        
                      
                      
                        ∑
                        
                          
                            c
                            j
                          
                          ⁢
                          
                            z
                            j
                            t
                          
                        
                      
                    
                  
                  )
                
                .
              
            
          
        
      
    
    
      
        (
        43
        )
      
    
  




f) Find a “small” D*s.t. zit+1(D*) that satisfies some “external” conditions (e.g., is a “positive” as a covariance matrix, one can use approximate representation Eq. (43) to find D*).

Gt({zit+1(D*)})≧Gt({zit+1(D)})  (44).

Remark: To find D* that makes a square covariance matrix positive one needs to solve a quadratic equation in D−1:





  
    
      
        
          
            
              
                (
                
                  
                    a
                    11
                  
                  +
                  
                    
                      D
                      
                        -
                        1
                      
                    
                    ⁢
                    
                      b
                      11
                    
                  
                
                )
              
              ⁢
              
                (
                
                  
                    a
                    22
                  
                  +
                  
                    
                      D
                      
                        -
                        1
                      
                    
                    ⁢
                    
                      b
                      22
                    
                  
                
                )
              
            
            -
            
              
                (
                
                  
                    a
                    12
                  
                  +
                  
                    
                      D
                      
                        -
                        1
                      
                    
                    ⁢
                    
                      b
                      12
                    
                  
                
                )
              
              2
            
          
          ≥
          0
        
        ⁢
        
          

        
        ⁢
        
          
            
              where
              ⁢
              
                
              
              ⁢
              
                a
                ij
              
            
            +
            
              
                D
                
                  -
                  1
                
              
              ⁢
              
                b
                ij
              
            
          
          =
          
            
              z
              i
              t
            
            +
            
              
                
                  ∑
                  
                    
                      c
                      i
                    
                    ⁢
                    
                      z
                      i
                      t
                    
                  
                
                D
              
              ⁢
              
                
                  (
                  
                    1
                    -
                    
                      
                        
                          z
                          i
                          t
                        
                        ⁢
                        
                          c
                          i
                        
                      
                      
                        ∑
                        
                          
                            c
                            j
                          
                          ⁢
                          
                            z
                            j
                            t
                          
                        
                      
                    
                  
                  )
                
                .
              
            
          
        
      
    
    
      
        (
        45
        )
      
    
  




g) xit+1=zit+1(D*)  (46).

Go to Step b).
For an even more general case, we consider the following unconstrained problem:

max{tilde over (Q)}(P,∥P∥1,λ)  (47),

where PεRL, λεR1. A recursion scheme for (47) includes the following:

    
    
        1. Find λt(P,∥P∥1) that satisfies
    
    






  
    
      
        
          
            ∂
            
              
                Q
                ~
              
              ⁡
              
                (
                
                  
                    P
                    t
                  
                  ,
                  
                    
                      
                      
                        P
                        t
                      
                      
                    
                    1
                  
                  ,
                  λ
                
                )
              
            
          
          
            ∂
            λ
          
        
        =
        0
      
    
    
      
        (
        48
        )
      
    
  




    
    
         as function of P,∥P∥1.
        2. Solve:

maxQ(P,∥P∥1)={tilde over (Q)}(Pt,∥Pt∥1,λt(P,∥P∥1))  (49)

    
    


Linear Problems: The benefit of using EBW for gradient computation is that can be used to solve optimization problems of large scale. For example, when trying to solve the following linear problem y=Hx+ε, LASSO was unable to solve this problem when the number of parameters x exceeded 100,000 whereas EBW was able to provide a solution. For example, FIG. 2 shows that EBW can reconstruct a noisy signal of 100,000 parameters. The squares in the figure are the true signal and the estimated signals are shown by crosses.
Non-Linear Problems: EBW can also be used to solve non-linear problems of the form y=ƒ(x)+ε. Here ƒ(x) is some non-linearity function over parameters x. One approach for solving this non-linear problem is by using an Extended Kalman Filter (EKF). FIGS. 3 and 4 respectfully show signal reconstruction using EKF and EBW. The squares are the true signal and the estimated signal are shown by lines. It can be seen that EBW does a much better job of signal reconstruction compared to EKF.
It should be understood that the optimization in accordance with the present principles provides a wide range of applications. In the example shown, signal reconstruction is performed using the optimization in accordance with the present embodiments. Signal reconstruction is applicable to many areas of endeavor, e.g., data classification, data compression, image processing, model probability density functions, speech recognition, telephony applications, financial market applications, etc. The signal reconstruction can involve a very large number of parameters in many such applications, e.g., in astronomy (when signals from whole space are received by a large number of antennas and a need to reconstruct special kinds of signals from noise generated by space radio signals exists), in fMRI analysis which can involve several hundred thousands of voxels, in telephone networks that can involve billions of Internet nodes where a need exists to analyze network traffic to detect patterns in certain situations (e.g., potential increase in terrorist activities), or in financial markets monitoring.
Referring to FIG. 5, a block diagram shows a system 100 for optimizing a function 110 in accordance with the present principles. System 100 is particularly useful for signal reconstruction applications. System 100 includes at least one processor 102 configured to determine a best gradient as a sparse direction in the function having a plurality of parameters. The sparse direction includes a direction that maximizes change of the function as determined by an optimization method where a sparse direction is defined from sparse regularized constraints. The processor 102 is configured to determine a best step size along the sparse direction by solving an optimization problem to find magnitudes of entries of direction that are constrained to the sparse direction. A memory storage device 104 includes memory storage media. The memory device 104 includes one or more programs 106. The program 106 is configured to recursively refine a stored solution for the function optimization.
The processor 102 and program 106 determine a size E of a directional ball which is employed as a constraint for determining the best direction in parameter space for the plurality of parameters. Alternatively, ε may be user selected. The processor 102 preferably performs an extended Baum Welch optimization 108 to identify the sparse direction. This is applied to the directional ball to find the sparse direction through a point of interest. It should be understood that other optimization methods may be employed. Given the sparse direction, the processor 102 and/or program also determines a larger sized ball R along the sparse direction and preferably performs the extended Baum-Welch method to identify a sub-optimal size inside R. Other optimizations methods may also be employed. The function 110 may include a function subject to a fractional normalized constrained inequality. The present principles are particularly useful if the function 110 includes at least 100,000 parameters.
Having described preferred embodiments of a system and method directional optimization via EBW (which are intended to be illustrative and not limiting), it is noted that modifications and variations can be made by persons skilled in the art in light of the above teachings. It is therefore to be understood that changes may be made in the particular embodiments disclosed which are within the scope of the invention as outlined by the appended claims. Having thus described aspects of the invention, with the details and particularity required by the patent laws, what is claimed and desired protected by Letters Patent is set forth in the appended claims.

## Claims

1. An optimization method, comprising:
determining a best gradient as a sparse direction in a function having a plurality of parameters, the sparse direction including a direction that maximizes change of the function as determined by an optimization method where a sparse direction is defined from sparse regularized constraints;
determining a best step size along the sparse direction by finding magnitudes of entries of direction by solving an optimization problem that is constrained to the sparse direction; and
recursively refining a solution to the optimization problem using a processor and storage media until a threshold constraint is met.

2. The method as recited in claim 1, wherein determining a best gradient includes identifying a size ε of a directional ball for determining the best direction for a directional optimization problem in parameter space of the plurality of parameters.

3. The method as recited in claim 2, wherein determining a best gradient includes given the sparse direction, determining a larger directional ball R along the sparse direction; and identifying a sub-optimal size inside R.

4. The method as recited in claim 1, wherein determining a best gradient includes performing an extended Baum Welch method to identify the sparse direction.

5. The method as recited in claim 4, wherein performing an extended Baum Welch (EBW) method to identify the sparse direction includes solving the EBW with the following:




  
    
      
        
          z
          i
          
            t
            +
            1
          
        
        =
        
          
            
              z
              i
              
                t
                +
                1
              
            
            ⁡
            
              (
              D
              )
            
          
          =
          
            
              
                (
                
                  
                    c
                    i
                  
                  +
                  D
                
                )
              
              ⁢
              
                z
                i
                t
              
            
            
              
                
                  ∑
                  j
                
                ⁢
                
                  
                    c
                    j
                  
                  ⁢
                  
                    z
                    j
                    t
                  
                
              
              +
              D
            
          
        
      
    
    
      
        (
        12
        )
      
    
  




where t, i and j are indexes,




  
    
      c
      j
    
    =
    
      
        
          c
          j
        
        ⁡
        
          (
          
            z
            t
          
          )
        
      
      =
      
        
          ∂
          
            G
            ⁡
            
              (
              
                z
                t
              
              )
            
          
        
        
          ∂
          
            z
            j
            t
          
        
      
    
  
  ,




D is a constant that controls growth of function values, G is a differentiable function in z, which converges to a local maximum of G(z) in the domain Σzi=1, zi≧0.

6. The method as recited in claim 1, wherein determining a best gradient includes computing the gradient as



  
    
      
        max
        v
      
      ⁢
      
        F
        ⁡
        
          (
          
            
              v
              0
            
            +
            v
          
          )
        
      
    
    -
    
      F
      ⁡
      
        (
        
          v
          0
        
        )
      
    
  




  
    
      s
      .
      t
      .
      
        
      
      ⁢
      
        
          
          
            v
            0
          
          
        
        1
      
    
    <
    ɛ
  



where F is the function, v is a direction, direction vector v belongs to a subspace that passes through the point v0, and a sparsity constraint includes ∥v∥1<ε where ε is a threshold for the sparsity constraint.

7. The method as recited in claim 1, wherein the function includes a function subject to a fractional normalized constrained inequality.

8. The method as recited in claim 7, further comprising:
adding a dumb non-negative variable to map the constrained inequality into a fractionally constrained multi-dimensional sphere;
transforming the fractional multi-dimensional sphere into a norm one sphere;
mapping the norm one sphere into a simplex of non-negative numbers; and
transforming an objective function over the norm one sphere into the objective function over the simplex.

9. A computer readable storage medium comprising a computer readable program for an optimization method, wherein the computer readable program when executed on a computer causes the computer to perform the steps of:
determining a best gradient as a sparse direction in a function having a plurality of parameters, the sparse direction including a direction that maximizes change of the function as determined by an optimization method where a sparse direction is defined from sparse regularized constraints;
determining a best step size along the sparse direction by finding magnitudes of entries of direction by solving an optimization problem that is constrained to the sparse direction; and
recursively refining a solution to the optimization problem until a threshold constraint is met.

10. The computer readable storage medium as recited in claim 9, wherein determining a best gradient includes identifying a size ε of a directional ball for determining the best direction for a directional optimization problem in parameter space of the plurality of parameters.

11. The computer readable storage medium as recited in claim 10, wherein determining a best gradient includes given the sparse direction, determining a larger directional ball R along the sparse direction; and identifying a sub-optimal size inside R.

12. The computer readable storage medium as recited in claim 9, wherein determining a best gradient includes performing an extended Baum Welch method to identify the sparse direction.

13. The computer readable storage medium as recited in claim 12, wherein performing an extended Baum Welch (EBW) method to identify the sparse direction includes solving the EBW with the following:




  
    
      
        
          z
          i
          
            t
            +
            1
          
        
        =
        
          
            
              z
              i
              
                t
                +
                1
              
            
            ⁡
            
              (
              D
              )
            
          
          =
          
            
              
                (
                
                  
                    c
                    i
                  
                  +
                  D
                
                )
              
              ⁢
              
                z
                i
                t
              
            
            
              
                
                  ∑
                  j
                
                ⁢
                
                  
                    c
                    j
                  
                  ⁢
                  
                    z
                    j
                    t
                  
                
              
              +
              D
            
          
        
      
    
    
      
        (
        12
        )
      
    
  




where t, i and j are indexes,




  
    
      c
      j
    
    =
    
      
        
          c
          j
        
        ⁡
        
          (
          
            z
            t
          
          )
        
      
      =
      
        
          ∂
          
            G
            ⁡
            
              (
              
                z
                t
              
              )
            
          
        
        
          ∂
          
            z
            j
            t
          
        
      
    
  
  ,




D is a constant that controls growth of function values, G is a differentiable function in z, which converges to a local maximum of G(z) in the domain Σzi=1, zi≦0.

14. The computer readable storage medium as recited in claim 9, wherein determining a best gradient includes computing the gradient as



  
    
      
        max
        v
      
      ⁢
      
        F
        ⁡
        
          (
          
            
              v
              0
            
            +
            v
          
          )
        
      
    
    -
    
      F
      ⁡
      
        (
        
          v
          0
        
        )
      
    
  




  
    
      s
      .
      t
      .
      
        
      
      ⁢
      
        
          
          
            v
            0
          
          
        
        1
      
    
    <
    ɛ
  



where F is the function, v is a direction, direction vector v belongs to a subspace that passes through the point v0, and a sparsity constraint includes ∥v∥1<ε where ε is a threshold for the sparsity constraint.

15. The computer readable storage medium as recited in claim 9, wherein the function includes a function subject to a fractional normalized constrained inequality.

16. The computer readable storage medium as recited in claim 15, further comprising:
adding a dumb non-negative variable to map the constrained inequality into a fractionally constrained multi-dimensional sphere;
transforming the fractional multi-dimensional sphere into a norm one sphere;
mapping the norm one sphere into a simplex of non-negative numbers; and
transforming an objective function over the norm one sphere into the objective function over the simplex.

17. A system for optimization of a function with a large number of parameters, comprising:
a processor configured to determine a best gradient as a sparse direction in a function having a plurality of parameters, the sparse direction including a direction that maximizes change of the function as determined by an optimization method where a sparse direction is defined from sparse regularized constraints;
the processor being configured to determine a best step size along the sparse direction by solving an optimization problem to find magnitudes of entries of direction that are constrained to the sparse direction; and
a memory storage device is configured to store a recursively refined optimization solution of the function.

18. The system as recited in claim 17, further comprising a directional ball of a size ε which is employed as a constraint for determining the best direction in parameter space for the plurality of parameters.

19. The system as recited in claim 18, wherein given the sparse direction, the processor determines a larger sized ball R along the sparse direction; and identifies a sub-optimal size inside R.

20. The system as recited in claim 17, wherein the processor performs an extended Baum Welch to identify the sparse direction.

21. The system as recited in claim 17, wherein a best gradient is determined by:



  
    
      
        max
        v
      
      ⁢
      
        F
        ⁡
        
          (
          
            
              v
              0
            
            +
            v
          
          )
        
      
    
    -
    
      F
      ⁡
      
        (
        
          v
          0
        
        )
      
    
  




  
    
      s
      .
      t
      .
      
        
      
      ⁢
      
        
          
          v
          
        
        1
      
    
    <
    ɛ
  



where F is the function, v is a direction, direction vector v belongs to a subspace that passes through the point v0 and a sparsity constraint includes ∥v∥1<ε where ε is a threshold for the sparsity constraint.

22. The system as recited in claim 21, wherein an extended Baum Welch (EBW) method identifies the sparse direction by solving the EBW with the following:




  
    
      
        
          z
          i
          
            t
            +
            1
          
        
        =
        
          
            
              z
              i
              
                t
                +
                1
              
            
            ⁡
            
              (
              D
              )
            
          
          =
          
            
              
                (
                
                  
                    c
                    i
                  
                  +
                  D
                
                )
              
              ⁢
              
                z
                i
                t
              
            
            
              
                
                  ∑
                  j
                
                ⁢
                
                  
                    c
                    j
                  
                  ⁢
                  
                    z
                    j
                    t
                  
                
              
              ⁢
              
                
              
              +
              D
            
          
        
      
    
    
      
        (
        12
        )
      
    
  




where t, i and j are indexes,




  
    
      c
      j
    
    =
    
      
        
          c
          j
        
        ⁡
        
          (
          
            z
            t
          
          )
        
      
      =
      
        
          ∂
          
            G
            ⁡
            
              (
              
                z
                t
              
              )
            
          
        
        
          ∂
          
            z
            j
            t
          
        
      
    
  
  ,




D is a constant that controls growth of function values, G is a differentiable function in z, which converges to a local maximum of G(z) in the domain Σzi=1, zi≧0.

23. The system as recited in claim 17, wherein the function includes a function subject to a fractional normalized constrained inequality.

24. The system as recited in claim 17, wherein the function includes at least 100,000 parameters.